Prev: HELP PLEASE!!
Next: remove black bar outline in bar3
From: monalisa on 11 May 2010 16:51
Sorry for late reply..thanks a ton !
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hs78l9$769$1(a)fred.mathworks.com>...
> 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
From: Bruno Luong on 11 May 2010 17:01
"monalisa " <mde1506(a)gmail.com> wrote in message <hscfvq$dmh$1(a)fred.mathworks.com>...
> Sorry for late reply..thanks a ton !

Smile, good to know you are still reading and it helps. It takes me some effort to make the code working right (check up to n = 20).

Bruno
From: monalisa on 11 May 2010 17:22
I learnt from the discussions :)...I really appreciate your effort...it took me some time to understand the algorithm. For my case, n=46 and its taking some time...but I am trying it right now.
Thanks again.

"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hscgij$krq$1(a)fred.mathworks.com>...
> "monalisa " <mde1506(a)gmail.com> wrote in message <hscfvq$dmh$1(a)fred.mathworks.com>...
> > Sorry for late reply..thanks a ton !
>
> Smile, good to know you are still reading and it helps. It takes me some effort to make the code working right (check up to n = 20).
>
> Bruno
From: Bruno Luong on 11 May 2010 17:34
"monalisa " <mde1506(a)gmail.com> wrote in message <hschps$ae3$1(a)fred.mathworks.com>...
> I learnt from the discussions :)...I really appreciate your effort...it took me some time to understand the algorithm. For my case, n=46 and its taking some time...but I am trying it right now.
> Thanks again.

If you try such largen, the weak point of the algorithm is perhaps the line of computing the NULL of Aeq matrix. Actually I have replace this part by a new function where I just submitted on FEX:
http://www.mathworks.com/matlabcentral/fileexchange/27550

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.';
TMat = sparse(p(:), idx(:), 1, n2, n2);

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

I = speye(n);

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

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

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

% SPNULL: http://www.mathworks.com/matlabcentral/fileexchange/27550
Q = spnull(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 = TMat*kron(EE,I);

M12 = M1+M2;

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

lbd = lbdref;
nresold = Inf;

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

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

TEMat = speye(n2) + TMat;

% Newton iteration
for i=1:maxiter

% Matrix for DX*X' + X*DX'
M34 = TEMat*kron(X,I);

% Group all together
M = M12+M34;

% 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 % for-loop
end
end % for-loop

if nres>nresold
X = XOLD;
end

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

end % LMsolve

% Bruno
From: monalisa on 11 May 2010 17:43
The solution is converging , but many x(i,j) are negative in the final X matrix.
does it mean one of the constraint condition is not working?


"monalisa " <mde1506(a)gmail.com> wrote in message <hschps$ae3$1(a)fred.mathworks.com>...
> I learnt from the discussions :)...I really appreciate your effort...it took me some time to understand the algorithm. For my case, n=46 and its taking some time...but I am trying it right now.
> Thanks again.
>
> "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hscgij$krq$1(a)fred.mathworks.com>...
> > "monalisa " <mde1506(a)gmail.com> wrote in message <hscfvq$dmh$1(a)fred.mathworks.com>...
> > > Sorry for late reply..thanks a ton !
> >
> > Smile, good to know you are still reading and it helps. It takes me some effort to make the code working right (check up to n = 20).
> >
> > Bruno
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: HELP PLEASE!!
Next: remove black bar outline in bar3