Solving Matrix equation with contraints
in [Matlab]
|
Prev: HELP PLEASE!!
Next: remove black bar outline in bar3
From: Bruno Luong on 11 May 2010 17:50 "monalisa " <mde1506(a)gmail.com> wrote in message <hscj19$pf5$1(a)fred.mathworks.com>... > 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? Do you want to impose X them positive? May be the easiest option is to use one of the function (e.g., FMINCON) in optimization toolbox. May be John can guide you, I do not have the license for this toolbox. Bruno
From: Bruno Luong on 11 May 2010 17:57 Sorry I overlooked your first post and the condition X>=0 is indeed required. Another option is just project X to the positive set after updating X % ... X = XOLD + DX; % new codes X = max(X,0); DX = X-XOLD; % .... Bruno
From: monalisa on 11 May 2010 20:12 I really appreciate your patience. When I am projecting X on the positive set,after updating in the final iteration, its giving me for some of the row sums > 1. "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hscjrg$g98$1(a)fred.mathworks.com>... > Sorry I overlooked your first post and the condition X>=0 is indeed required. Another option is just project X to the positive set after updating X > > % ... > X = XOLD + DX; > % new codes > X = max(X,0); > DX = X-XOLD; > % > ... > > Bruno
From: Bruno Luong on 12 May 2010 01:24 "monalisa " <mde1506(a)gmail.com> wrote in message <hscrom$pvd$1(a)fred.mathworks.com>... > I really appreciate your patience. When I am projecting X on the positive set,after updating in the final iteration, its giving me for some of the row sums > 1. Sorry, I forget to take into account other constraints. What I can offer is making the following change: find the Newton step-size (called rho) so that the condition X>=0 is prevailed. I make two alterations in the code, finding the best maximum step-size and in the test accordingly for breaking the Newton's for-loop. Remember that the more you increase the size, the more the problem becomes unstable. Indeed, the matrix equation looks like a non-linear first kind Fredholm's integral equation, which is a prototype of ill-posed problem. Any numerical solution will have a limitation. Here is the new code. 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 % X >= 0 % % 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-10; % Levenberg-Marquardt parameters lbdref = 0.01; v = 2; % 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); TEMat = speye(n2) + TMat; % Project on the affine space X(:) = X0(:) + P*(X(:)-X0(:)); % 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]); % Find the step maximum step size so as % to prevent X to be < 0 after updating DX(1:n+1:n) = 0; rho = -XOLD./DX; rho(DX>=0) = Inf; rho = min(min(rho(:)),1); % Update X X = XOLD + rho*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 = rho*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 12 May 2010 13:15
Thanks a ton....I am unfamiliar with Fredholm's integral equation, but I will read it up. Thanks again. "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <hsde1k$jd9$1(a)fred.mathworks.com>... > "monalisa " <mde1506(a)gmail.com> wrote in message <hscrom$pvd$1(a)fred.mathworks.com>... > > I really appreciate your patience. When I am projecting X on the positive set,after updating in the final iteration, its giving me for some of the row sums > 1. > > Sorry, I forget to take into account other constraints. What I can offer is making the following change: find the Newton step-size (called rho) so that the condition X>=0 is prevailed. I make two alterations in the code, finding the best maximum step-size and in the test accordingly for breaking the Newton's for-loop. > > Remember that the more you increase the size, the more the problem becomes unstable. Indeed, the matrix equation looks like a non-linear first kind Fredholm's integral equation, which is a prototype of ill-posed problem. Any numerical solution will have a limitation. > > Here is the new code. > > 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 > % X >= 0 > % > % 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-10; > % Levenberg-Marquardt parameters > lbdref = 0.01; > v = 2; % 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); > > TEMat = speye(n2) + TMat; > > % Project on the affine space > X(:) = X0(:) + P*(X(:)-X0(:)); > > % 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]); > > % Find the step maximum step size so as > % to prevent X to be < 0 after updating > DX(1:n+1:n) = 0; > rho = -XOLD./DX; > rho(DX>=0) = Inf; > rho = min(min(rho(:)),1); > % Update X > X = XOLD + rho*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 = rho*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 |