11

151行151列の行列がありAます。これは相関行列であるため1、主対角線上に s があり、主対角線の上下に値が繰り返されます。各行/列は人を表します。

与えられた整数に対して、要素の総和を最小化する相関行列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最適な行列は、人物 1 と 4 を含む2行2列の恒等行列であることは明らかですが、反復スキームによれば、反復の最初のフェーズで人物 1 を追い出してしまうため、アルゴリズムは次のような解を作成します。最適ではありません。私は常に最適な解を生成するプログラムを書きました。n または k が小さい場合はうまく機能しますが、151 行 x 151 列から最適な75行75 列の行列を作成しようとすると、151行列 私のプログラムが終了するには数十億年かかることに気付きました。

これらのn を選択する kの問題は、再計算を回避する動的プログラミングのアプローチで解決できる場合があることを漠然と思い出しましたが、これを解決する方法がわかりません。

他に選択肢がない場合、または最良のプログラムでも正確な解を生成するのに 1 週​​間以上かかる場合は、速度のために精度を犠牲にしてもかまいません。ただし、正確な解が生成される場合は、プログラムを最大 1 週間実行させていただきます。

プログラムが合理的な時間枠内で行列を最適化できない場合、この特定の種類のn を選択した kタスクを合理的な時間枠内で解決できない理由を説明する答えを受け入れます。

4x4 相関行列

4

4 に答える 4

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

は、2 つの親を結合する新しい個体を生成する機能です。

この時点で、問題を定義できます。

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 に答える
3

近似解を見つけるにはいくつかのアプローチがあります (例: 緩和された問題に対する二次計画法や貪欲な探索) が、正確な解を見つけることは NP 困難な問題です。

免責事項: 私は 2 項二次計画法の専門家ではありません。より高度なアルゴリズムについては学術文献を参照してください。

数学的に同等の定式化:

あなたの問題は次と同等です：

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) と呼ばれます。バイナリ バージョンは、Binary Quadratic Programming (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は元の問題の解の下限であることに注意してください。

あなたの正確な問題を解決するためのソフトウェア

このリンクをチェックして、混合整数二次計画法 (MIQP) のセクションに移動してください。Gurobi はあなたのタイプの問題を解決できるように見えます。別のソルバーのリストはこちらです。

于 2015-11-21T20:01:20.973 に答える
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 に答える