11

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

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

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

相关矩阵

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

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

我明白那个sum(A(:))将汇总矩阵中的所有内容。但是，我非常不清楚如何搜索最小可能的答案。

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

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

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

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

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

4x4 相关矩阵

4

4 回答 4

3

如果我理解你的问题陈述，你有一个N x N矩阵M（它恰好是一个相关矩阵），你希望找到整数n其中 2 <= n < N，一个n x n矩阵m最小化我用f (m)表示的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 的索引）

于 2015-11-27T16:53:51.807 回答
3

There are several approaches to finding an approximate solution (eg. quadratic programming on relaxed problem or greedy search), but finding the exact solution is an NP-hard problem.

Disclaimer: I'm not an expert on binary quadratic programming, and you may want to consult the academic literature for more sophisticated algorithms.

Mathematically equivalent formulation:

Your problem is equivalent to:

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

This is a quadratic programming problem where the vector x is restricted to taking only binary values. Quadratic programming where the domain is restricted to a set of discrete values is called mixed integer quadratic programming (MIQP). The binary version is sometimes called Binary Quadratic Programming (BQP). The last restriction, that x is binary, makes the problem substantially more difficult; it destroys the problem's convexity!

Quick and dirty approach to finding an approximate answer:

If you don't need a precise solution, something to play around with might be a relaxed version of the problem: drop the binary constraint. If you drop the constraint that x(i) is either 1 or 0 for all i, then the problem becomes a trivial convex optimization problem and can be solved nearly instantaneously (eg. by Matlab's quadprog). You could try removing entries that, on the relaxed problem, quadprog assigns the lowest values in the x vector, but this does not truly solve the original problem!

Note also that the relaxed problem gives you a lower bound on the optimal value of the original problem. If your discretized version of the solution to the relaxed problem leads to a value for the objective function close to the lower bound, there may be a sense in which this ad-hoc solution can't be that far off from the true solution.

To solve the relaxed problem, you might try something like:

% 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;

I'm curious how good x_approx is? This doesn't solve your problem, but it might not be horrible! Note that f_relax is a lower bound on the solution to the original problem.

Software to solve your exact problem

You should 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.

于 2015-11-21T20:01:20.973 回答
3

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

我从你的测试用例开始：

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)只有所需人员的矩阵也是如此。

于 2015-11-25T17:02:58.283 回答
1

根据 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
于 2015-11-28T08:22:48.333 回答