MATLAB：找到最小化矩阵元素总和的矩阵的缩写版本
MATLAB: Find abbreviated version of matrix that minimises sum of matrix elements
我有一个 151-by-151 矩阵 A。它是一个相关矩阵，因此主对角线上有 1，主对角线上下有重复值。每个row/column代表一个人。
对于给定的整数 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 的最佳矩阵是 2-by-2 涉及第 1 个人和第 4 个人的单位矩阵，但根据迭代方案我会在迭代的第一阶段踢出 Person 1，因此该算法会得出一个不是最优的解决方案。我写了一个总是生成最优解的程序，当 n 或 k 很小时它运行良好，但是当试图做出最优 75-by-75[=来自 151-by-151 矩阵的 64=] 矩阵 我意识到我的程序需要数十亿年才能终止。
我依稀记得有时候这些 n choose k 问题可以用避免重新计算的动态规划方法来解决，但我不知道如何解决这个问题，并且谷歌搜索也没有启发我。
如果没有其他选择，我愿意牺牲精度来换取速度，否则最好的程序将需要一周以上的时间才能生成精确的解决方案。但是，如果程序能够生成精确的解决方案，我很乐意让它 运行 长达一周。
如果程序不可能在合理的时间范围内优化矩阵，那么我会接受一个答案，解释为什么 n 选择 k 这种特定类型的任务不能'无法在合理的时间内解决。
有几种方法可以找到近似解（例如，松弛问题的二次规划或贪婪搜索），但 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 向量中分配最低值的条目，但这 并没有 真正解决原始问题！
另请注意，松弛问题为您提供了原始问题最优值的下限。如果您对松弛问题的解决方案的离散版本导致 objective 函数的值接近下限，则可能存在这种临时解决方案不能远离下限的感觉真正的解决方案。
要解决放松的问题，您可以尝试这样的方法：
% 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 and go down to the section on Mixed Integer Quadratic Programming (MIQP). It looks to me that Gurobi can solve problems of your type. Another list of solvers is here.
这是使用遗传算法的近似解。
我从你的测试用例开始：
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 菜单 select Import Problem... 并在打开的 window 中选择 gaproblem2。
现在，运行 工具并等待迭代停止。
gatool 使您能够更改数百个参数，因此您可以在 selected 输出中以速度换取精度。
生成的向量是您必须保留在原始矩阵中的索引列表，因此 A(garesults.x,garesults.x) 是仅包含所需人员的矩阵。
如果我理解你的问题陈述，你有一个 N x N 矩阵 M（它恰好是一个相关矩阵），你希望找到整数 n where 2 <= n < N, 一个 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)
根据 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
我有一个 151-by-151 矩阵 A。它是一个相关矩阵，因此主对角线上有 1，主对角线上下有重复值。每个row/column代表一个人。
对于给定的整数 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 的最佳矩阵是 2-by-2 涉及第 1 个人和第 4 个人的单位矩阵，但根据迭代方案我会在迭代的第一阶段踢出 Person 1，因此该算法会得出一个不是最优的解决方案。我写了一个总是生成最优解的程序，当 n 或 k 很小时它运行良好，但是当试图做出最优 75-by-75[=来自 151-by-151 矩阵的 64=] 矩阵 我意识到我的程序需要数十亿年才能终止。
我依稀记得有时候这些 n choose k 问题可以用避免重新计算的动态规划方法来解决，但我不知道如何解决这个问题，并且谷歌搜索也没有启发我。
如果没有其他选择，我愿意牺牲精度来换取速度，否则最好的程序将需要一周以上的时间才能生成精确的解决方案。但是，如果程序能够生成精确的解决方案，我很乐意让它 运行 长达一周。
如果程序不可能在合理的时间范围内优化矩阵，那么我会接受一个答案，解释为什么 n 选择 k 这种特定类型的任务不能'无法在合理的时间内解决。
有几种方法可以找到近似解（例如，松弛问题的二次规划或贪婪搜索），但 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 向量中分配最低值的条目，但这 并没有 真正解决原始问题！
另请注意，松弛问题为您提供了原始问题最优值的下限。如果您对松弛问题的解决方案的离散版本导致 objective 函数的值接近下限，则可能存在这种临时解决方案不能远离下限的感觉真正的解决方案。
要解决放松的问题，您可以尝试这样的方法：
% 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 and go down to the section on Mixed Integer Quadratic Programming (MIQP). It looks to me that Gurobi can solve problems of your type. Another list of solvers is here.
这是使用遗传算法的近似解。
我从你的测试用例开始：
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 菜单 select Import Problem... 并在打开的 window 中选择 gaproblem2。
现在，运行 工具并等待迭代停止。
gatool 使您能够更改数百个参数，因此您可以在 selected 输出中以速度换取精度。
生成的向量是您必须保留在原始矩阵中的索引列表，因此 A(garesults.x,garesults.x) 是仅包含所需人员的矩阵。
如果我理解你的问题陈述，你有一个 N x N 矩阵 M（它恰好是一个相关矩阵），你希望找到整数 n where 2 <= n < N, 一个 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)
根据 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