mchol - % if (j >= 2), if (j < n),

Info iconThis preview shows pages 1–3. Sign up to view the full content.

View Full Document Right Arrow Icon
% % [L,D,E,pneg]=mchol1(G) % % Given a symmetric matrix G, find a matrix E of "small" norm and c % L, and D such that G+E is Positive Definite, and % % G+E = L*D*L' % % Also, calculate a direction pneg, such that if G is not PD, then % % pneg'*G*pneg < 0 % % Note that if G is PD, then the routine will return pneg=[]. % % Reference: Gill, Murray, and Wright, "Practical Optimization", p111. % Author: Brian Borchers (borchers@nmt.edu) % function [L,D,E,pneg]=mchol(G) % % n gives the size of the matrix. % n=size(G,1); % % gamma, zi, nu, and beta2 are quantities used by the algorithm. % gamma=max(diag(G)); zi=max(max(G-diag(diag(G)))); nu=max([1,sqrt(n^2-1)]); beta2=max([gamma, zi/nu, 1.0E-15]); % % Initialize diag(C) to diag(G). % C=diag(diag(G)); % % Loop through, calculating column j of L for j=1:n % L=zeros(n); D=zeros(n); E=zeros(n); for j=1:n, bb=[1:j-1]; ee=[j+1:n]; % % Calculate the jth row of L. % if (j > 1), L(j,bb)=C(j,bb)./diag(D(bb,bb))'; end; % % Update the jth column of C.
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
Background image of page 3
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: % if (j &gt;= 2), if (j &lt; n), C(ee,j)=G(ee,j)-(L(j,bb)*C(ee,bb)')'; end; else C(ee,j)=G(ee,j); end; % % Update theta. % if (j == n) theta(j)=0; else theta(j)=max(abs(C(ee,j))); end; % % Update D % D(j,j)=max([eps,abs(C(j,j)),theta(j)^2/beta2]'); % % Update E. % E(j,j)=D(j,j)-C(j,j); % % Update C again. .. % %%%%%%%% M.Zibulevsky: begin of changes, old version is commented %%%%%%%%%%%%% %for i=j+1:n, % C(i,i)=C(i,i)-C(i,j)^2/D(j,j); %end; ind=[j*(n+1)+1 : n+1 : n*n]'; C(ind)=C(ind)-(1/D(j,j))*C(ee,j).^2; end; % % Put 1's on the diagonal of L % %for j=1:n, % L(j,j)=1; %end; ind=[1 : n+1 : n*n]'; L(ind)=1; %%%%%%%%%%%%%%%%%%%%%%%% M.Zibulevsky: end of changes %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % if needed, find a descent direction. % if ((nargout == 4) &amp; (min(diag(C)) &lt; 0.0)) [m,col]=min(diag(C)); rhs=zeros(n,1); rhs(col)=1; pneg=L'\rhs; else pneg=; end; return...
View Full Document

Page1 / 3

mchol - % if (j &amp;amp;gt;= 2), if (j &amp;amp;lt; n),

This preview shows document pages 1 - 3. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online