MATLAB：查找最小化矩阵元素总和的矩阵的缩写版本答案

【问题标题】：MATLAB: Find abbreviated version of matrix that minimises sum of matrix elementsMATLAB：查找最小化矩阵元素总和的矩阵的缩写版本
【发布时间】：2016-02-17 18:09:39
【问题描述】：

我有一个 151-by-151 矩阵 A。它是一个相关矩阵，因此主对角线上有1s，主对角线上方和下方有重复值。每行/列代表一个人。

对于给定的整数n，我将通过将人们踢出去来减少矩阵的大小，这样我就剩下一个n-by-n 相关矩阵，它可以最小化元素的总和。除了获得缩写矩阵，我还需要知道应该从原始矩阵中引导出来的人的行号（或他们的列号 - 他们将是相同的数字）。

我以A = tril(A) 为起点，它将从相关矩阵中删除多余的非对角元素。

所以，如果n = 4 并且我们有上面假设的 5-by-5 矩阵，那么很明显应该将人 5 踢出矩阵，因为那个人贡献了很多非常高的相关性。

很明显，第 1 个人不应该被踢出去，因为那个人贡献了很多负相关，从而降低了矩阵元素的总和。

我知道sum(A(:)) 将对矩阵中的所有内容求和。但是，我非常不清楚如何搜索尽可能少的答案。

我注意到一个类似的问题Finding sub-matrix with minimum elementwise sum，它有一个蛮力解决方案作为公认的答案。虽然该答案在那里工作正常，但对于 151-by-151 矩阵来说是不切实际的。

编辑：我曾想过迭代，但我认为这并不能真正最小化缩减矩阵中的元素总和。下面我有一个 4-by-4 粗体相关矩阵，边上有行和列的总和。很明显，n = 2 的最佳矩阵是涉及人员 1 和 4 的 2-by-2 单位矩阵，但根据迭代方案，我会踢出人 1 处于迭代的第一阶段，因此算法得出的解决方案不是最优的。我编写了一个程序，它总是生成最优解，当 n 或 k 很小时它运行良好，但是当试图从151-by-151 矩阵我意识到我的程序需要数十亿年才能终止。

我隐约记得，有时这些 n 选择 k 问题可以通过避免重新计算事物的动态编程方法来解决，但我不知道如何解决这个问题，谷歌搜索也没有启发我.

如果没有其他选择，我愿意牺牲精度来换取速度，否则最好的程序需要一周多的时间才能生成精确的解决方案。不过，如果程序能够生成精确的解决方案，我很乐意让程序运行长达一周。

如果程序无法在合理的时间范围内优化矩阵，那么我会接受一个答案，解释为什么这种特定类型的 n 选择 k 任务无法在合理的时间内解决时间范围。

【问题讨论】：

sum(A, 2) 返回每一行的总和。
您是否能够扩展该代码如何帮助找到 n-by-n 矩阵的缩写版本，该矩阵的总和最小化矩阵元素？
让我看看我是否准确地理解了这个问题。问题是选择一个向量x 以最小化x'*S*x，其中S 是给定的对称正定矩阵，x 受制于x 的条目是1 或@ 987654342@ 和x 的元素总和为n。对吗？
This link 可能会有所帮助。另外，您是否愿意为了速度而牺牲精度？
我编辑了 OP 以阐明精度与速度的问题。您对问题的理解是正确的，除了 S 是相关矩阵，因此必然是半正定的，不一定是正定的。

标签： matlab optimization quadratic-programming

【解决方案1】：

有几种方法可以找到近似解（例如，放松问题或贪婪搜索的二次规划），但finding the exact solution is an NP-hard problem。

免责声明：我不是二元二次规划方面的专家，您可能需要查阅学术文献以了解更复杂的算法。

数学等价公式：

你的问题相当于：

For some symmetric, positive semi-definite matrix S

minimize (over vector x)  x'*S*x

             subject to     0 <= x(i) <= 1        for all i
                            sum(x)==n
                            x(i) is either 1 or 0 for all i

这是一个二次规划问题，其中向量x 仅限于采用二进制值。将域限制为一组离散值的二次规划称为混合整数二次规划 (MIQP)。二进制版本有时称为二进制二次规划 (BQP)。最后一个限制，即x 是二进制的，使得问题实质上变得更加困难；它破坏了问题的凸性！

找到近似答案的快速而肮脏的方法：

如果您不需要精确的解决方案，可以尝试解决问题的一个轻松版本：删除二元约束。如果你放弃x(i) is either 1 or 0 for all i 的约束，那么问题就变成了一个微不足道的凸优化问题，几乎可以立即解决（例如，通过Matlab 的quadprog）。您可以尝试删除在宽松问题上 quadprog 分配 x 向量中最低值的条目，但这并没有真正解决原始问题！

另请注意，松弛问题为您提供了原始问题最优值的下限。如果您对松弛问题的离散化版本的解决方案导致目标函数的值接近下限，那么在某种意义上，这种临时解决方案可能与真正的解决方案相差不远。

要解决轻松的问题，您可以尝试以下方法：

% k is number of observations to drop
n = size(S, 1);
Aeq = ones(1,n)
beq = n-k;
[x_relax, f_relax] = quadprog(S, zeros(n, 1), [], [], Aeq, beq, zeros(n, 1), ones(n, 1));
f_relax = f_relax * 2;  % Quadprog solves .5 * x' * S * x... so mult by 2
temp = sort(x_relax);
cutoff = temp(k);
x_approx = ones(n, 1);
x_approx(x_relax <= cutoff) = 0;
f_approx = x_approx' * S * x_approx;

我很好奇 x_approx 有多好？这不能解决你的问题，但它可能并不可怕！请注意，f_relax 是原始问题解决方案的下限。

解决您确切问题的软件

您应该check out this link 并转到混合整数二次规划 (MIQP) 部分。在我看来，Gurobi 可以解决您类型的问题。另一个list of solvers is here。

【讨论】：

确实问题和你刚才说的一样，除了一件事。相关矩阵是半正定的而不是正定的。例如，参见stats.stackexchange.com/questions/182875/…。这会改变事情吗？
@user1205901 基本相同。如果随机变量的某些线性组合与另一个完全相关，则可以得到零特征值（因此是半定而不是定）。
0 <= x(i) <= 1 for all i 有错字吗？ S 中的元素必须介于 0 和 1 之间，而 x 中的元素必须是 0 或 1。
如果你有最后一个，显然第一个是多余的。我只是把它放在那里，所以我可以干净利落地争辩说后一个问题只是前者的限制较少的版本。
签出我的链接去 Gurobi 了吗？这可能是你想要的。你说解决方案是针对pastebin中的那个？我认为 [1;1;1;0;0] 是您在 patebin 中发布的 S 矩阵的解决方案。

【解决方案2】：

这是使用遗传算法的近似解。

我从你的测试用例开始：

data_points = 10; % How many data points will be generated for each person, in order to create the correlation matrix.
num_people = 25; % Number of people initially.
to_keep = 13; % Number of people to be kept in the correlation matrix.
to_drop = num_people - to_keep; % Number of people to drop from the correlation matrix.
num_comparisons = 100; % Number of times to compare the iterative and optimization techniques.
for j = 1:data_points
    rand_dat(j,:) = 1 + 2.*randn(num_people,1); % Generate random data.
end
A = corr(rand_dat);

然后我定义了进化遗传算法所需的函数：

function individuals = user1205901individuals(nvars, FitnessFcn, gaoptions, num_people)

individuals = zeros(num_people,gaoptions.PopulationSize);
for cnt=1:gaoptions.PopulationSize
    individuals(:,cnt)=randperm(num_people);
end

individuals = individuals(1:nvars,:)';

是个体生成函数。

function fitness = user1205901fitness(ind, A)

fitness = sum(sum(A(ind,ind)));

是适应度评价函数

function offspring = user1205901mutations(parents, options, nvars, FitnessFcn, state, thisScore, thisPopulation, num_people)

offspring=zeros(length(parents),nvars);
for cnt=1:length(parents)
    original = thisPopulation(parents(cnt),:);
    extraneus = setdiff(1:num_people, original);
    original(fix(rand()*nvars)+1) = extraneus(fix(rand()*(num_people-nvars))+1);
    offspring(cnt,:)=original;
end

是改变个体的功能

function children = user1205901crossover(parents, options, nvars, FitnessFcn, unused, thisPopulation)

children=zeros(length(parents)/2,nvars);
cnt = 1;
for cnt1=1:2:length(parents)
    cnt2=cnt1+1;
        male = thisPopulation(parents(cnt1),:);
        female = thisPopulation(parents(cnt2),:);
        child = union(male, female);
        child = child(randperm(length(child)));
        child = child(1:nvars);
        children(cnt,:)=child;
        cnt = cnt + 1;

end

是生成耦合两个父母的新个体的函数。

此时您可以定义您的问题：

gaproblem2.fitnessfcn=@(idx)user1205901fitness(idx,A)
gaproblem2.nvars = to_keep
gaproblem2.options = gaoptions()
gaproblem2.options.PopulationSize=40
gaproblem2.options.EliteCount=10
gaproblem2.options.CrossoverFraction=0.1
gaproblem2.options.StallGenLimit=inf
gaproblem2.options.CreationFcn= @(nvars,FitnessFcn,gaoptions)user1205901individuals(nvars,FitnessFcn,gaoptions,num_people)
gaproblem2.options.CrossoverFcn= @(parents,options,nvars,FitnessFcn,unused,thisPopulation)user1205901crossover(parents,options,nvars,FitnessFcn,unused,thisPopulation)
gaproblem2.options.MutationFcn=@(parents, options, nvars, FitnessFcn, state, thisScore, thisPopulation) user1205901mutations(parents, options, nvars, FitnessFcn, state, thisScore, thisPopulation, num_people)
gaproblem2.options.Vectorized='off'

打开遗传算法工具

gatool

从File 菜单中选择Import Problem... 并在打开的窗口中选择gaproblem2。

现在，运行该工具并等待迭代停止。

gatool 可让您更改数百个参数，因此您可以用速度换取所选输出的精度。

生成的向量是您必须保留在原始矩阵中的索引列表，因此A(garesults.x,garesults.x) 是仅包含所需人员的矩阵。

【讨论】：

【解决方案3】：

如果我理解你的问题陈述，你有一个 N x N 矩阵 M （恰好是一个相关矩阵），并且您希望找到整数 n 其中 2 n N, a n x n 矩阵 m 最小化 m 的所有元素的总和，我表示为 f(m)?

在 Matlab 中，获取矩阵的子矩阵相当容易和快速（例如参见 Removing rows and columns from matrix in Matlab），并且函数 f 相对便宜评估 n = 151。那么为什么你不能在下面的程序中实现一个算法来动态解决这个问题，我已经勾勒出伪代码：

function reduceM(M, n){

    m = M

    for (ii = N to n+1) {

        for (jj = 1 to ii) {

            val(jj) = f(m) where mhas column and row jj removed, f(X) being summation over all elements of X

        }

        JJ(ii) = jj s.t. val(jj) is smallest

        m = m updated by removing column and row JJ(ii)   
    }
}

最后你会得到一个维度为 n 的 m，它是你的问题的解决方案和一个向量 JJ，其中包含在每次迭代中删除的索引（你应该能够轻松地将这些转换回适用于完整的索引矩阵 M)

【讨论】：

这个答案完全误解了问题的复杂性。从这个问题中，“我注意到一个类似的问题'找到具有最小元素和的子矩阵'，它有一个强力解决方案作为公认的答案。虽然该答案在那里工作得很好，但对于 151×151 矩阵来说是不切实际的。”如果用户想要保留 140 行/列（消除 11），则需要 nchoosek(151, 140) = 大约 10^29 个不同的矩阵来检查！那是行不通的。 ...solves this backwards dynamically... 含糊不清。这个问题是一个非常难的问题，老实说这个答案没有帮助。
感谢您的评论，因为您看到的复杂性来自您选择的解决方案方法，将其视为一个组合问题。
我们的不同之处在于我建议我们通过反向归纳方法来解决它。给定一个大小为 M 的矩阵，如果我们只是希望将维度减少 1，即 n = M -1，我们可以很容易地找到问题的解决方案。如果我们将这个过程应用 M - n 次，我们会得到一个对 n 最优的答案，并且我们可以记录沿途移除的“人”。
不，那不一定是你给出的答案！！想象两个问题：（1）将维度减少 1 或（2）将维度减少 2。它们可能涉及删除完全不同的行/列。
如果您理解我的建议，您会意识到我们并没有对所有可能的解决方案矩阵进行 10^29 次检查，而更像是大约 x = 151+150+149...n 次，因为在每个迭代你必须检查尽可能多的可能解决方案作为起始矩阵的维度

【解决方案4】：

根据 Matthew Gunn 的建议以及 Gurobi 论坛上的一些建议，我提出了以下功能。它似乎工作得很好。

我会给它答案，但如果有人能提出更好的代码，我会从这个答案中删除勾号并将其放在他们的答案上。

function [ values ] = the_optimal_method( CM , num_to_keep)
%the_iterative_method Takes correlation matrix CM and number to keep, returns list of people who should be kicked out 

N = size(CM,1);  

clear model;  
names = strseq('x',[1:N]);  
model.varnames = names;  
model.Q = sparse(CM); % Gurobi needs a sparse matrix as input  
model.A = sparse(ones(1,N));  
model.obj = zeros(1,N);  
model.rhs = num_to_keep;  
model.sense = '=';  
model.vtype = 'B';

gurobi_write(model, 'qp.mps');

results = gurobi(model);

values = results.x;

end

【讨论】：