Prev: HELP PLEASE!!
Next: remove black bar outline in bar3
From: Bruno Luong on 9 May 2010 08:30
"John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <hs61va$k2g$1(a)fred.mathworks.com>...
> "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!!!!!

Sorry John, what did I replace? Let me allow to quote OP

["I am trying to solve AX=B;"]

Then Wikipedia page: http://en.wikipedia.org/wiki/Equation_solving

[ In mathematics, to solve an equation is to find what values (numbers, functions, sets, etc.) fulfill a condition stated in the form of an equation (two expressions related by equality). These expressions contain one or more unknowns, which are free variables for which values are sought that cause the condition to be fulfilled. To be precise, what is sought are often not necessarily actual values, but, more in general, mathematical expressions. A solution of the equation is an assignment of expressions to the unknowns that satisfies the equation; in other words, expressions such that, when they are substituted for the unknowns, the equation becomes a tautology (a provably true statement).
...

It is also possible that the task is to find a solution, among possibly many, that is best in some respect. Problems of that nature are called optimization problems; solving an optimization problem is generally *not* referred to as "equation solving".
]

The last sentence is what in my mind. OP asked to *solve* an equation (not to approximate it). I just put three equations he stated in matrix form. I did *not* replace anything.

Bruno
From: John D'Errico on 9 May 2010 10:26
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hs69sh$d2p$1(a)fred.mathworks.com>...
> "John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <hs61va$k2g$1(a)fred.mathworks.com>...
> > "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!!!!!
>
> Sorry John, what did I replace? Let me allow to quote OP
>
> ["I am trying to solve AX=B;"]

SUBJECT to a set of equality constraints.

When you append the equality constraints to the linear
system A*X = B, the system becomes overdetermined.
It no longer has an exact solution in general.

Therefore, the equality constraints will no longer be
exactly satisfied. The problems are very different, although
one might argue that the difference is a subtle one, the
difference is important.

Surely you understand that there is a difference here?

For example, solve the linear system A*x = b.

A = rand(2,3)
b = rand(2,1)

A =
0.814723686393179 0.126986816293506 0.63235924622541
0.905791937075619 0.913375856139019 0.0975404049994095
b =
0.278498218867048
0.546881519204984

A\b
ans =
0.293942721840442
0.307245445468768
0

This is an underdetemined problem, as you well know.

However, suppose we choose to solve that linear
system subject to a set of EQUALITY constraints?

Aeq = rand(2,3);
beq = rand(2,1);

Aeq
Aeq =
0.957506835434298 0.157613081677548 0.957166948242946
0.964888535199277 0.970592781760616 0.485375648722841

beq
beq =
0.8002804688888
0.141886338627215

Here we solve the linear system A*x = b, subject to
the exact satisfaction of the linear equality constraints.

x0 = lsqlin(A,b,[],[],Aeq,beq)
x0 =
1.45935801310315
-1.08173016982432
-0.445658909531652

What does this do in terms of the original system?

A*x0 - b
ans =
0.491293318876076
-0.256502768131462

Aeq*x0 - beq
ans =
0
1.38777878078145e-16

See that the equality constraints are exactly that -
equality constraints. They are satisfied exactly. See
what happens when we combine the two matrices
into one and just use backslash.

x1 = [A;Aeq]\[b;beq]
x1 =
1.07802458175086
-0.58604426595662
-0.394212388430651

See that the equality constraints are not satisfied.
In fact, they are no longer equality constraints.

Aeq*x1 - beq
ans =
-0.237759874646762
0.138135792868382

A*x1 - b
ans =
0.276090198057855
-0.144145864262541

John
From: Bruno Luong on 9 May 2010 11:01
"John D'Errico" <woodchips(a)rochester.rr.com> wrote in message <hs6glr$mf1$1(a)fred.mathworks.com>...

>
> A*x0 - b
> ans =
> 0.491293318876076
> -0.256502768131462
>
> Aeq*x0 - beq
> ans =
> 0
> 1.38777878078145e-16
>
> See that the equality constraints are exactly that -
> equality constraints. They are satisfied exactly. See
> what happens when we combine the two matrices
> into one and just use backslash.
>

But John, your example does NOT provide the solution for A*x0 = b, regardless which methods (since (A*x0 - b) is NOT zeros).

OP asked to SOLVE A*x = b, i.e., A*x - b = 0 at the numerical precision.

Bruno
From: Bruno Luong on 9 May 2010 12:09
Here is a more robust code using levenberg-marquardt damping:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% script LMsolve.m
% Solve for X the system
% E*(X+X')+C*(X*X') = D
% diag(X) = 0
% sum(X,2) = 1
%
% All matrices are square and real and have the same size (n x n)

%% Data
n=4;

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

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

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

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

%% Engine
n = length(D);
n2 = n*n;
idx = zeros(n);
idx(:) = 1:numel(idx);

% Newton pamameters
maxiter = 500;
tol = 1e-6;
% levenberg-marquardt parameters
lbdref = 0.01;
v = 2;

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 orthogonal span <rows Aeq>

% New matrix equation: EE*(X+X')+(X*X') = DD
EE = C\E;
DD = C\D;

% Matrix for EE*DX
M1 = kron(speye(n),EE);

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

converge = false;
resfun = @(X) DD - (EE*(X+X') + X*X');

lbd = lbdref;
nresold = Inf;

% Random starting solution
X = zeros(n);

% Newton iteration
for i=1:maxiter

% Project on the affine space
X(:) = X0(:) + P*(X(:)-X0(:));

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

% Matrix for DX*X'
M4 = kron(X',speye(n))';

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

% Diagonal of Levenberg-marquardt
LM = sqrt(sum(M.^2,1));
LM = spdiags(LM(:), 0, n2, n2);

% Rhs
B = resfun(X);
RHS = [B(:); zeros(n2,1)];

% Keep track old matrix
XOLD = X;

% LM iteration on damping factor lbd
while true
J = [M; lbd*LM];

DX = J \ RHS;
DX = P*DX;
DX = reshape(DX,[n n]);

% Update X
X(:) = XOLD(:) + DX(:);

residual = resfun(X);
nres = norm(residual(:));
if nres<nresold || lbd>1e3
nresold = nres;
if lbd>1e-3
lbd = lbd/v;
end
break % while-loop
else
lbd = lbd*v;
end
end % LM inner-loop
nDX = norm(DX(:));
if nDX<tol
converge = true;
break
end
end % for-loop

X = XOLD;

%% Check
if converge
fprintf('#iterations = %d\n', i);
else
warning('not converge');
end
fprintf('X output:\n')
disp(X)
fprintf('E*(X+X'') + C*X*X''\n')
disp(E*(X+X') + C*(X*X'))

%% Test on command line

>> LMsolve
X input:
0 0.2284 0.4930 0.2786
0.7168 0 0.1256 0.1576
0.1499 0.6263 0 0.2238
0.4942 0.4217 0.0841 0

D
-0.1010 -1.1638 -1.0334 -0.8023
0.5220 0.0579 -1.1870 -0.9846
-0.2692 -2.2210 -0.4314 -0.9152
-0.2210 -0.7911 -1.0095 -1.0191

#iterations = 10
X output:
0 0.2284 0.4930 0.2786
0.7168 0 0.1256 0.1576
0.1499 0.6263 0 0.2238
0.4942 0.4217 0.0841 0

E*(X+X') + C*X*X'
-0.1010 -1.1638 -1.0334 -0.8023
0.5220 0.0579 -1.1870 -0.9846
-0.2692 -2.2210 -0.4314 -0.9152
-0.2210 -0.7911 -1.0095 -1.0191

% Bruno
From: Bruno Luong on 9 May 2010 17:15
I make an important change in a way to perform the internal projection. The new code is below:

function X = LMsolve(E, C, D)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X = LMsolve(E, C, D)
% Solve for X the system
% E*(X+X')+C*(X*X') = D
% diag(X) = 0
% sum(X,2) = 1
%
% All matrices are square and real and have the same size (n x n)
% Levenberg-Marquardt method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Engine, try to find X from E, C, D
n = length(D);
n2 = n*n;
idx = reshape(1:n2,[n n]);
p = idx.';
transm = sparse(p(:), idx(:), 1, n2, n2);

% Newton pamameters
maxiter = 1000;
tol = 1e-6;
% Levenberg-Marquardt parameters
lbdref = 0.01;
v = 4;

I = speye(n);

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

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

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

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

% Reformulate matrix equation: EE*(X+X')+(X*X') = DD
% and other constraints
EE = C\E;
DD = C\D;

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

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

M12 = M1+M2;

converge = false;
resfun = @(X) DD - (EE*(X+X') + X*X');

lbd = lbdref;
nresold = Inf;

% "Random" starting solution
X = zeros(n);

% Newton iteration
for i=1:maxiter

% Project on the affine space
X(:) = X0(:) + P*(X(:)-X0(:));

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

% Matrix for X*DX'
M3 = transm*M4;

% Group all together
M = M12+M3+M4;

% We look for DX := Q*DY
% Such that M*DX = RHS
M = M*Q;

% Diagonal of Levenberg-marquardt damping
LM = sqrt(sum(M.^2,1));
m = size(M,2);
LM = sparse(1:m, 1:m, LM, m, m);

% Rhs
B = resfun(X);
RHS = [B(:); zeros(m,1)];

% Keep track old matrix
XOLD = X;

% LM iteration on damping factor lbd
while true
J = [M; lbd*LM];

DY = J\RHS;
DX = Q*DY;
DX = reshape(DX,[n n]);

% Update X
X(:) = XOLD(:) + DX(:);

residual = resfun(X);
nres = norm(residual(:));
if nres<=nresold || lbd>1e3
nresold = nres;
if lbd>1e-6
lbd = lbd/v;
end
break % while-loop
else
lbd = lbd*v;
end
end % LM inner-loop
nDX = norm(DX(:));
if nDX<tol
converge = true;
break
end
end % for-loop

X = XOLD;

%% Check
if converge
fprintf('#iterations = %d\n', i);
else
warning('not yet converge');
end

end % LMsolve

%%%%%

% Test LMSolve, a Levenberg-Marquardt solver
%% Data
n=4;

E = rand(n)+diag(rand(1,n));
C = rand(n)+diag(rand(1,n));

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

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

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

%% Command line
X = LMsolve(E, C, D);

fprintf('------------------------\n')
fprintf('Levenberg-Marquardt: X output:\n')
disp(X)
fprintf('E*(X+X'') + C*X*X'':\n');
DX = E*(X+X') + C*(X*X');
disp(DX);

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