Solves AX + XA = C when A = A', as a pseudo-inverse, given eig(A). function X = lyapunov_symmetric_eig(V, lambda, C) function X = lyapunov_symmetric_eig(V, lambda, C, tol) Same as lyapunov_symmetric(A, C, [tol]), where A is symmetric, its eigenvalue decomposition [V, lambda] = eig(A, 'vector') is provided as input directly, and C is a single matrix of the same size as A. See also: lyapunov_symmetric sylvester lyap sylvester_nocheck
0001 function X = lyapunov_symmetric_eig(V, lambda, C, tol) 0002 % Solves AX + XA = C when A = A', as a pseudo-inverse, given eig(A). 0003 % 0004 % function X = lyapunov_symmetric_eig(V, lambda, C) 0005 % function X = lyapunov_symmetric_eig(V, lambda, C, tol) 0006 % 0007 % Same as lyapunov_symmetric(A, C, [tol]), where A is symmetric, its 0008 % eigenvalue decomposition [V, lambda] = eig(A, 'vector') is provided as 0009 % input directly, and C is a single matrix of the same size as A. 0010 % 0011 % See also: lyapunov_symmetric sylvester lyap sylvester_nocheck 0012 0013 % This file is part of Manopt: www.manopt.org. 0014 % Original author: Nicolas Boumal, Aug. 31, 2018. 0015 % Contributors: 0016 % Change log: 0017 0018 % AX + XA = C is equivalent to DY + YD = M with 0019 % Y = V'XV, M = V'CV and D = diag(lambda). 0020 M = V'*C*V; 0021 0022 % W(i, j) = lambda(i) + lambda(j) 0023 W = bsxfun(@plus, lambda, lambda'); 0024 0025 % Normally, the solution Y is simply this: 0026 Y = M ./ W; 0027 0028 % But this may involve divisions by (almost) 0 in certain places. 0029 % Thus, we go for a pseudo-inverse. 0030 absW = abs(W); 0031 if ~exist('tol', 'var') || isempty(tol) 0032 tol = numel(C)*eps(max(absW(:))); % similar to pinv tolerance 0033 end 0034 Y(absW <= tol) = 0; 0035 0036 % Undo the change of variable 0037 X = V*Y*V'; 0038 0039 end