Поскольку M является двоичным, мы можем думать об этом как о проблеме графа. i, j в {1..n} соответствуют узлам, а M(i,j) указывает, существует ли неориентированный край, соединяющий их.
С A, B, C, D являются непересекающимися, что немного упрощает проблему. Мы можем подойти к решению проблемы поэтапно:
- Найти все
(c,d), для которых существует a такое, что M(a,c) < M(a,d). Назовем этот набор CD_lt_a (подмножество C * D такое, что для некоторого выполняется неравенство «меньше».
- Найти все
(c,d), для которых существует a такой, что M(a,c) <= M(a,d), и позвоните по этому набору CD_le_a.
- Повторите для
b, образуя CD_lt_b для M(b,d) < M(b,c) и CD_le_b для M(b,d)<=M(b,c).
- Один из способов удовлетворить общее неравенство -
M(a,c) < M(a,d) и M(b,d) <= M(b,c), поэтому мы можем посмотреть на пересечение CD_lt_a и CD_le_b.
- Другой способ - если
M(a,c) <= M(a,d) и M(b,d) < M(b,c), поэтому перейдите на перекресток CD_le_a и CD_lt_b.
- С
(c,d) известно, мы можем вернуться назад и найти (a,b).
И вот моя реализация:
% 0. Some preliminaries
% Get the size of each set
mA = numel(A); mB = numel(B); mC = numel(C); mD = numel(D);
% 1. Find all (c,d) for which there exists a such that M(a,c) < M(a,d)
CA_linked = M(C,A);
AD_linked = M(A,D);
CA_not_linked = ~CA_linked;
% Multiplying these matrices tells us, for each (c,d), how many nodes
% in A satisfy this M(a,c)<M(a,d) inequality
% Ugh, we need to cast to double to use the matrix multiplication
CD_lt_a = (CA_not_linked * double(AD_linked)) > 0;
% 2. For M(a,c) <= M(a,d), check that the converse is false for some a
AD_not_linked = ~AD_linked;
CD_le_a = (CA_linked * double(AD_not_linked)) < mA;
% 3. Repeat for b
CB_linked = M(C,B);
BD_linked = M(B,D);
CD_lt_b = (CB_linked * double(~BD_linked)) > 0;
CD_le_b = (~CB_linked * double(BD_linked)) < mB;
% 4. Find the intersection of CD_lt_a and CD_le_b - this is one way
% to satisfy the inequality M(a,c)+M(b,d) < M(a,d)+M(b,c)
CD_satisfy_ineq_1 = CD_lt_a & CD_le_b;
% 5. The other way to satisfy the inequality is CD_le_a & CD_lt_b
CD_satisfy_ineq_2 = CD_le_a & CD_lt_b;
inequality_feasible = any(CD_satisfy_ineq_1(:) | CD_satisfy_ineq_2(:));
Обратите внимание, что вы можете остановиться здесь, если осуществимость ваша единственная забота. Сложность составляет A*C*D + B*C*D, что лучше, чем худший случай A*B*C*D сложность цикла for. Однако раннее завершение означает, что ваши вложенные для циклов все еще могут быть быстрее в определенных случаях.
В следующем блоке кода перечислены все a,b,c,d, которые удовлетворяют неравенству. Он не очень хорошо оптимизирован (он присоединяется к матрице изнутри цикла), поэтому он может быть довольно медленным, если есть много результатов.
% 6. With (c,d) known, find a and b
% We can define these functions to help us search
find_a_lt = @(c,d) find(CA_not_linked(c,:)' & AD_linked(:,d));
find_a_le = @(c,d) find(CA_not_linked(c,:)' | AD_linked(:,d));
find_b_lt = @(c,d) find(CB_linked(c,:)' & ~BD_linked(:,d));
find_b_le = @(c,d) find(CB_linked(c,:)' | ~BD_linked(:,d));
% I'm gonna assume there aren't too many results, so I will be appending
% to an array inside of a for loop. Bad for performance, but maybe a bit
% more readable for a StackOverflow answer.
results = zeros(0,4);
% Find those that satisfy it the first way
[c_list,d_list] = find(CD_satisfy_ineq_1);
for ii = 1:numel(c_list)
c = c_list(ii); d = d_list(ii);
a = find_a_lt(c,d);
b = find_b_le(c,d);
% a,b might be vectors, in which case all combos are valid
% Many ways to find all combos, gonna use ndgrid()
[a,b] = ndgrid(a,b);
% Append these to the growing list of results
abcd = [a(:), b(:), repmat([c d],[numel(a),1])];
results = [results; abcd];
end
% Repeat for the second way
[c_list,d_list] = find(CD_satisfy_ineq_2);
for ii = 1:numel(c_list)
c = c_list(ii); d = d_list(ii);
a = find_a_le(c,d);
b = find_b_lt(c,d);
% a,b might be vectors, in which case all combos are valid
% Many ways to find all combos, gonna use ndgrid()
[a,b] = ndgrid(a,b);
% Append these to the growing list of results
abcd = [a(:), b(:), repmat([c d],[numel(a),1])];
results = [results; abcd];
end
% Remove duplicates
results = unique(results, 'rows');
% And actually these a,b,c,d will be indices into A,B,C,D because they
% were obtained from calling find() on submatrices of M.
if ~isempty(results)
results(:,1) = A(results(:,1));
results(:,2) = B(results(:,2));
results(:,3) = C(results(:,3));
results(:,4) = D(results(:,4));
end
Я испытал это на следующий тест:
m = 1000;
A = (1:m); B = A(end)+(1:m); C = B(end)+(1:m); D = C(end)+(1:m);
M = rand(D(end),D(end)) < 1e-6; M = M | M';
Мне нравится думать, что первая часть (см если неравенство является допустимым для любого а, Ь, с, d) работал очень хорошо , Другие векторизованные ответы (которые используют ndgrid или combvec для перечисления всех комбинаций a, b, c, d) потребуют 8 терабайт памяти для проблемы такого размера!
Но я бы не рекомендовал запускать вторую часть (перечисляя все результаты), когда есть несколько сотен c,d, которые удовлетворяют неравенству, потому что это будет довольно медленно.
P.S.Я знаю, что уже ответил, но этот ответ касался векторизации таких циклов в целом и менее специфичен для вашей конкретной проблемы.
P.P.S. Этот вид напоминает мне об stable marriage problem. Возможно, некоторые из этих ссылок будут содержать алгоритмы, относящиеся к вашей проблеме. Я подозреваю, что истинный алгоритм, основанный на графике, может, вероятно, достичь наихудшей сложности, так как это дополнительно предлагает раннее завершение. Но я думаю, что было бы сложно эффективно реализовать алгоритм на основе графов в MATLAB.
P.P.P.S. Если вы хотите только одно из возможных решений, вы можете упростить шаг 6, чтобы вернуть только одно значение, например.
find_a_lt = @(c,d) find(CA_not_linked(c,:)' & AD_linked(:,d), 1, 'first');
find_a_le = @(c,d) find(CA_not_linked(c,:)' | AD_linked(:,d), 1, 'first');
find_b_lt = @(c,d) find(CB_linked(c,:)' & ~BD_linked(:,d), 1, 'first');
find_b_le = @(c,d) find(CB_linked(c,:)' | ~BD_linked(:,d), 1, 'first');
if any(CD_satisfy_ineq_1)
[c,d] = find(CD_satisfy_ineq_1, 1, 'first');
a = find_a_lt(c,d);
b = find_a_le(c,d);
result = [A(a), B(b), C(c), D(d)];
elseif any(CD_satisfy_ineq_2)
[c,d] = find(CD_satisfy_ineq_2, 1, 'first');
a = find_a_le(c,d);
b = find_a_lt(c,d);
result = [A(a), B(b), C(c), D(d)];
else
result = zeros(0,4);
end
Вы можете использовать ответы от этого вопроса, чтобы получить комбинированную матрицу http://stackoverflow.com/questions/4165859/generate-all-possible-combinations-of-the-elements-of-some-vectors -cartesian-pr –
Пожалуйста, определите свои массивы в коде Matlab, а не на словах. –