MATLAB: найдите сокращенную версию матрицы, которая минимизирует сумму матричных элементов

Ответ 1

Работая над предложением Мэтью Ганна, а также некоторыми советами на форумах 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

Ответ 2

Существует несколько подходов к поиску приближенного решения (например, квадратичное программирование по расслабленной проблеме или жадному поиску), но найти точное решение - это NP-hard проблема.

Отказ от ответственности: я не специалист по двоичному квадратичному программированию, и вы можете обратиться к академической литературе для более сложных алгоритмов.

Математически эквивалентная формулировка:

Ваша проблема эквивалентна:

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, но это действительно не решает исходную проблему!

Заметим также, что релаксированная проблема дает вам нижнюю границу оптимального значения исходной задачи. Если ваша дискретизированная версия решения ослабленной проблемы приводит к значению целевой функции, близкой к нижней границе, может возникнуть смысл, в котором это ad-hoc решение не может быть настолько далеким от истинного решения.

Чтобы решить расслабленную проблему, вы можете попробовать что-то вроде:

% 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). Мне кажется, что Гуроби может решить проблемы вашего типа. Другой список решателей здесь.

Ответ 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

- это функция генерации новой индивидуальной связи двух родителей.

На этом этапе вы можете определить свою проблему:

[email protected](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)
[email protected](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) - это матрица с нужными лицами.

Ответ 4

Если я понял, что вы задаете проблему, у вас есть матрица N x N M (которая оказывается корреляционной матрицей) и вы хотите найти для integer n, где 2 <= n N, n x n матрица m, которая минимизирует сумму по всем элементам m, который я обозначаю f (m)?

В Matlab довольно легко и быстро получить подматрицу матрицы (см., например, Удаление строк и столбцов из матрицы в 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)   
    }
}

В итоге вы получаете m размера n, который является решением вашей проблемы, и вектор JJ, который содержит индексы, удаленные на каждой итерации (вы должны легко преобразовать их обратно в индексы, применимые к полной матрица M)