Solving Matrix equation with contraints [Matlab]

Prev: HELP PLEASE!!
Next: remove black bar outline in bar3

From: monalisa sen on 8 May 2010 11:33

"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hs3k1k$76b$1(a)fred.mathworks.com>...
> Actually since all the equation are equalities, no simple linear inversion (no need fancy toolbox) should do:
>
> % John's example
> n = 4;
> A = [ 16 2 3 13;
> 5 11 10 8;
> 9 7 6 12;
> 4 14 15 1];
> B = [4.78082343328984 9.22416092385747 11.3133835756658 8.68163206718693; ...
> 8.71780822242756 5.02919066456483 8.44394288910161 11.809058223906; ...
> 7.34746997932943 7.70521084234125 9.64129536853806 9.30602380979126; ...
> 8.89183816258422 1.19610039052821 7.7213261373564 16.1907353095312];
>
> Aeq = [kron(speye(size(A)),A);
> sparse(1:n,1:n+1:n^2,1,n,n*n);
> kron(ones(1,n),speye(n))];
>
> beq = [B(:);
> zeros(n,1);
> ones(n,1)];
>
> X = reshape(Aeq\beq,[n n]);
>
> disp(X)
>
> % Bruno

Thanks a lot John and Bruno. I really appreciate your help. Both my A and B matrix are full rank, but B can be real or complex matrix. How would I go about if it is complex.

thanks

From: monalisa sen on 8 May 2010 11:44

if I have to solve for X with the same constraints as above, but now the equation is
E(X+X')+C(XX')=D; all these matrices E,C,D are full rank and real.

any help/ suggestion would be great.

From: Bruno Luong on 8 May 2010 17:24

From: John D'Errico on 9 May 2010 06:15

"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hs3k1k$76b$1(a)fred.mathworks.com>...
> Actually since all the equation are equalities, no simple linear inversion (no need fancy toolbox) should do:
>

NO! I strongly disagree with this statement.

There is a very significant difference between what
Bruno shows and what I showed how to do.

When you solve a linear system of equations subject to
equality constraints as did I, you force the equality
constraints to hold exactly true. Thus, solve the linear
system

A*x = y

subject to the linear equality constraints

Aeq*x = beq

When you do so with a tool designed to solve that
linear system, you force the constraints in Aeq to
be satisfied EXACTLY, while the equations in A*x = y
are satisfied as well as you can subject to the equality
constraints. You are then minimizing

norm(A*x - y)

What Bruno has suggested doing is to replace this
equality constrained problem with a very different
one!!!!!

He tells you to form the new linear system

[A;Aeq] * x = [y;beq]

and to solve that system as a single least squares
problem. The equality constraints are now lumped
in with the rest of A, allowing them to be only
approximately satisfied!!!!! You are now solving
the problem of minimizing

norm([A;Aeq]*x - [y;beq])

This is a VERY different thing. A very different
solution, and I'll claim not an appropriate one
given the problem statement. That it is an easier
one to solve in MATLAB since it does not require
lsqlin does not make it a better solution.

As long as you have an exact solution, as I did
when I made up data in my test case, it did not
matter. The two methods must yield the same
result. However, in the real world, when your data
is imperfect, there is a significant difference!!!!!

John

From: Bruno Luong on 9 May 2010 07:30

"monalisa sen" <msen1981(a)gmail.com> wrote in message <hs40s5$7df$1(a)fred.mathworks.com>...
> if I have to solve for X with the same constraints as above, but now the equation is
> E(X+X')+C(XX')=D; all these matrices E,C,D are full rank and real.
>
> any help/ suggestion would be great.

Despite what I wrote, I can't resist to make a code to solve of the equation:

% script newtonsolve.m
%% Data
n=4;

E = 10*randn(n);
C = 10*randn(n);

% Generate X that saitisfies constraints
X = rand(n);
X(1:n+1:end) = 0;
X = diag(1./sum(X,2))*X;

% rhs
D = C*X*X'+E*(X+X');

fprintf('X input:\n')
disp(X);
fprintf('D\n')
disp(D);

%% Engine
idx = zeros(n);
idx(:) = 1:numel(idx);

maxiter = 100;
tol = 1e-10;

Aeq = [sparse(1:n,1:n+1:n^2,1,n,n*n);
kron(ones(1,n),speye(n))];

beq = [zeros(n,1);
ones(n,1)];

X0 = reshape(Aeq \ beq, [n n]);

Q = null(full(Aeq));
P=Q*Q'; % Projection operator on span <Aeq>

EE = C\E;
DD = C\D;

% Matrix for EE*DX
M1 = kron(eye(size(EE)),EE);

% Matrix for EE*DX'
M2 = kron(EE',eye(size(EE)))';
M2 = M2(idx',:);

% Random starting solution
X = zeros(n);

lhsfun = @(X) (EE*(X+X') + X*X');
% Newton iteration
for i=1:maxiter

% Matrix for X*DX'
M3 = kron(X',eye(size(X)))';
M3 = M3(idx',:);

% Matrix for DX*X'
M4 = kron(X',eye(size(X')))';

% Group all together
M = M1+M2+M3+M4;

% Rhs
B = DD - lhsfun(X);

DX = reshape(pinv(M)*B(:),[n n]);

XOLD = X;
% Update X, and make sure it satisfies the constraints
X(:) = X0(:) + P*(X(:)+DX(:)-X0(:));
res = norm(XOLD-X,'fro');
if res<tol
break
end
end

%% Check
fprintf('#iterations = %d\n', i);
fprintf('X output:\n')
disp(X)
fprintf('E*(X+X'') + C*X*X''\n')
disp(E*(X+X') + C*X*X')

%%%%%%%%%%%%%%%%%%%%%

>> newtonsolve
X input:
0 0.200011057410700 0.435143156125102 0.364845786464198
0.041699642636209 0 0.838590892356224 0.119709465007567
0.632003764772581 0.038824206504272 0 0.329172028723146
0.283411262822000 0.400519151217284 0.316069585960716 0

D
8.437875454824908 8.231567016110517 25.378926447605572 20.767614997066133
-32.576607481949722 -23.167874664568505 -23.733560823288471 -17.317195377377494
-21.109645575796083 -19.144084118959430 -28.679022904527024 -21.219115873571241
27.951333799118892 23.738997187380996 30.296712983976700 35.870650389416355

#iterations = 10
X output:
0 0.200011057410682 0.435143156125108 0.364845786464209
0.041699642636213 0 0.838590892356238 0.119709465007549
0.632003764772585 0.038824206504273 0 0.329172028723141
0.283411262821972 0.400519151217300 0.316069585960729 0

E*(X+X') + C*X*X'
8.437875454824923 8.231567016110521 25.378926447605522 20.767614997065859
-32.576607481949466 -23.167874664568380 -23.733560823288894 -17.317195377377303
-21.109645575795909 -19.144084118959562 -28.679022904527226 -21.219115873571123
27.951333799119027 23.738997187380903 30.296712983977219 35.870650389416028

% Bruno

First | Prev | Next | Last
Pages: 1 2 3 4 5
Prev: HELP PLEASE!!
Next: remove black bar outline in bar3