MATLAB：查找缩写版本的矩阵，使矩阵元素之和最小化

Question

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

给定整数 n，我将寻求通过剔除一些人来减小矩阵的大小，从而留下一个最小化元素总和的 n-by-n 相关系数矩阵。除了获得缩写矩阵外，我还需要知道应该从原始矩阵中剔除哪些人的行号（或它们的列号 - 它们将是同一个编号）。

作为起点，我取 A = tril(A)，这将从相关系数矩阵中删除冗余的非对角线元素。

因此，如果 n = 4，并且我们拥有上述假设的 5-by-5 矩阵，则非常明确的是应该从矩阵中剔除第五个人，因为该人贡献了很多非常高的相关性。

另外，显然不应该剔除第一个人，因为该人贡献了很多负的相关性，从而降低了矩阵元素的总和。

我知道 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 的矩阵生成一个最优的 75-by-75 矩阵时，我意识到我的程序将需要数十亿年才能终止。

我模糊地记得有时候这些 n choose k 问题可以通过动态规划方法来解决，避免重新计算，但我无法解决这个问题，谷歌也没有给我启示。

如果没有其他选择，我愿意为速度而牺牲精度，或者最好的程序将需要超过一周的时间来生成精确的解决方案。然而，如果一个程序能够生成精确的解决方案，我可以让它运行多达一周的时间。

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

- user1205901 - Слава Україні

---

sum(A, 2) 返回每行的总和。 - sco1

你能否详细说明一下这段代码如何帮助找到最小化矩阵元素和的n乘以n缩写版本的矩阵？ - user1205901 - Слава Україні

让我看看是否准确理解了问题。问题是选择一个向量 x，使得 x'*S*x 最小化，其中 S 是给定的对称正定矩阵，而 x 受到以下约束：x 的元素只能为 1 或 0，且 x 的元素之和为 n。是这样吗？ - Matthew Gunn

这个链接可能有帮助：http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-972-algebraic-techniques-and-semidefinite-optimization-spring-2006/lecture-notes/lecture_03.pdf。另外，你愿意为速度而牺牲精确性吗？ - Matthew Gunn

我编辑了原帖以澄清精度与速度问题。你对问题的理解是正确的，除了S是一个相关矩阵，因此必须是半正定的，而不一定是正定的。 - user1205901 - Слава Україні

Answer 1

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

我从您的测试案例开始：

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

“generate_offspring()”函数用于通过将两个父代进行交叉配对来生成新的个体。

在这一步，您可以定义您的问题：

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

从 文件 菜单中选择 导入问题... ，然后在打开的窗口中选择 gaproblem2。

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

gatool 可以让你改变数百个参数，因此你可以在所选输出中权衡速度和精度。

结果向量是您必须在原始矩阵中保留的索引列表，这样 A(garesults.x,garesults.x) 就是只包含所需人员的矩阵。

- NicolaSysnet

Answer 2

如果我理解你的问题陈述正确的话，你有一个 N x N 的矩阵 M（这恰好是一个相关矩阵），并且你希望找到整数n（其中2 <= n < N），一个n x n 的矩阵 m，使得所有元素的总和最小，我将其表示为f(m)？

在Matlab中，获得矩阵的子矩阵非常容易和快速（例如请参见Removing rows and columns from matrix in Matlab），对于 n = 151 来说，函数f的评估相对便宜。那么，为什么你不能实现下面用伪代码草绘出来的程序来动态地反向解决这个问题的算法呢？

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的索引）。

- Jimmy

---

感谢您的评论，您认为问题的复杂性来自于您选择的解决方法，将其视为组合问题。 - Jimmy

写你的代码。你的贪心算法不一定会收敛于全局最优解。句号。 - Matthew Gunn

我认错了。Matthew Gunn是正确的，我的解决方案是错误的。请忽略。 - Jimmy

没关系。这个问题表面看起来很简单... 直到我开始阅读一些内容才意识到它出奇的困难！ - Matthew Gunn

这就是我喜欢数学的原因 :) 我仍然认为有一种方法可以设计出一个高效的解决方案，希望有人能想出来！ - Jimmy

Answer 3

有几种方法可以找到近似解（例如，在放松问题上进行二次规划或贪婪搜索），但是找到确切的解是一个NP难问题。

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

数学等价形式：

您的问题等价于：

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)对于所有i都是1或0的约束条件，则问题变成一个微不足道的凸优化问题，可以通过几乎瞬间解决（例如，使用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是原始问题解的下界。

---

软件解决您的确切问题

您应该查看此链接并转到混合整数二次规划（MIQP）部分。在我看来，Gurobi可以解决您类型的问题。另一个求解器列表在这里。

- Matthew Gunn

---

这确实是一个问题，除了一点。相关矩阵是半正定而不是正定的。例如，请参见http://stats.stackexchange.com/questions/182875/is-every-correlation-matrix-positive-definite。这会改变事情吗？ - user1205901 - Слава Україні

@user1205901 基本上是同样的情况。如果一些随机变量的线性组合与另一个变量完全相关，则可能会得到零特征值（因此是半正定而不是正定）。 - Matthew Gunn

“0 <= x(i) <= 1 for all i” 中有错别字吗？实际上，S 中的元素必须介于 0 和 1 之间，而 x 中的元素必须是 0 或 1。 - user1205901 - Слава Україні

如果你有最后一个，第一个显然是多余的。我只是把它放在那里，这样我可以干净地证明后一个问题只是前一个问题的一个不太受限制的版本。 - Matthew Gunn

你看过我的Gurobi链接了吗？这可能是你想要的。那么你认为pastebin中的那个解决方案是什么？我认为[1;1;1;0;0]是你在pastebin上发布的S矩阵的解决方案。 - Matthew Gunn

Answer 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

- user1205901 - Слава Україні