Задача из теории графов

[Alex Krylov]

А можно вот так (для поиска вещественных решений, используем Matlab и парсер Yalmip):

Используется синтаксис Matlab M

x=sdpvar(46,1);
constr=[x>=0,A*x==s];
sol=optimize(constr, x(end));
value(x)

%Пример вещественного решения, такого, что sum(x(1:N))=x(N+1)=19
%x = 1.5193    0.7122    0.7122    0.7122    0.7122    1.0232    0.3427    0.3427    0.4254    0.4254    0.3427
%      0.3427    0.4254    0.4254    0.2927    0.3427    0.3427    0.4254    0.4254    0.2927    0.2927    0.3427
%      0.3427    0.4254    0.4254    0.2927    0.2927    0.2927    0.3427    0.3427    0.4254    0.4254    0.2927
%      0.2927    0.2927    0.2927    0.3427    0.3427    0.4254    0.4254    0.2927    0.2927    0.2927    0.2927
%      0.2927   19.0000
 

А затем можно вот так (для поиска двоичного решения):

Используется синтаксис Matlab M

x=binvar(45,1)
constr=[A(1:10,1:N)*x==s(1:10)];
sol=optimize(constr, sum(x));
value(x)

%Пример двоичного решения, дающего в сумме 19:
%x = 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1;
%sum(x)=19; 
 

[Stensen]

Спасибо, разбираюсь

[Alex Krylov]

Если обратиться к самому 1-му куску кода Matlab, то смысл там такой.

Пусть $x_{ij}, i,j=1..10$ - матрица смежности нашего графа. Эта матрица симметрична и имеет нулевую диагональ. Т.е. в этой матрице $(10^2-10)/2$ уникальных элементов. Пусть это элементы нижне-треугольной части матрицы, т.е. $x_{ij}, i,j=1..10, i>j$. Эту нижнетреугольную матрицу можно векторизовать и ввести сквознную линейную нумерацию: $\begin{bmatrix} *&* &* \\ x_{21}&* &* \\ x_{31}&x_{32} &* \end{bmatrix}\to\begin{bmatrix} x_{21}& x_{31} &x_{32}\end{bmatrix}$.
Тогда элементу $(i,j)$ матрицы смежности соответсвует уникальный элемент с линейным индексом $ind$:

Используется синтаксис Matlab M

        
q = max(i,j);
k = min(i,j);
ind = (q^2-3*q)/2+k+1;

Т.е. "x" - в вышеприведенном коде Matlab это (построчно) векторизованная нижнетреугольная часть матрицы смежности с дополнительным добавочным (46-м) элементом в конце. Этот добавочный 46-й элемент должен быть равен сумме предыдущих 45 элементов:

Используется синтаксис Matlab M

x(46)=sum(x(1:45));

Далее в виде СЛАУ $A x = s$ формализуется система линейных ограничений, приведенных в условии задачи. Последняя строка матрицы A (имеющая вид [-ones(1,45), 1]) и последний нулевой элемент вектора s используются для формализации вышеприведенного условия:

Используется синтаксис Matlab M

x(46)=sum(x(1:45));

Что из себя представляет общее решение неоднородной СЛАУ $A x=s$?
А вот что:

Используется синтаксис Matlab M

[U,S,V]=svd(A);

Spinv = S.';
for i=1:size(S,1) 
    Spinv(i,i)=1/Spinv(i,i);
end

pinvA = V*Spinv*U.';% псевдообратная матрица Мура-Пенроуза
sol=pinvA*s;
sol(end)
%Общее решение x неоднородной СЛАУ A*x=s имеет вид
% x = pinvA*s + null(A)*c;
% где null(A) - базис нуль-пространства матрицы A
% c - вектор-столбец произвольных коэффициентов.
 

Замечаем, что у всех векторов, образующих базис нуль-пространства, последний элемент равен НУЛЮ! Т.е. вклад нуль-пространства в общее решение x(end) равен нулю. Т.е. сумма уникальных элементов матрицы смежности определяется только частным решением СЛАУ и не зависит от ядра A.

[Alex Krylov]

Можно кстати говоря эту задачу решать через исключение по Гауссу, т.е. через LU-разложение. Это фактически символьный метод.

Используется синтаксис Matlab M

[L,U]=lu(A(1:10,1:45));
%U(:,1:35)*x(1:35) + U(:,36:45)*x(36:45) = L^-1*s(1:10)
%U(:,36:45)*x(36:45)=L^-1*s(1:10) - U(:,1:35)*x(1:35)
%x(36:45) = U(:,36:45)^-1*(L^-1*s(1:10) - U(:,1:35)*x(1:35))
%x(36:45) = U(:,36:45)^-1*L^-1*s(1:10) - U(:,36:45)^-1*U(:,1:35)*x(1:35)

%sum_x(1:35) = ones(1,35)*x(1:35)
%sum_x(36:45) = ones(1,10)*U(:,36:45)^-1*L^-1*s(1:10) - ones(1,10)*U(:,36:45)^-1*U(:,1:35)*x(1:35)
%sum_x = sum_x(1:35)+sum_x(36:45) = ones(1,10)*U(:,36:45)^-1*L^-1*s(1:10) +
% + (ones(1,35)-ones(1,10)*U(:,36:45)^-1*U(:,1:35))*x(1:35)
%Т.к. ones(1,35)-ones(1,10)*U(:,36:45)^-1*U(:,1:35) == 0,
%то sum_x = ones(1,10)*U(:,36:45)^-1*L^-1*s(1:10)

ones(1,10)*U(:,36:45)^-1*L^-1*s(1:10) % равно 19
ones(1,35)-ones(1,10)*U(:,36:45)^-1*U(:,1:35) % вектор-строка из нулей