-
Notifications
You must be signed in to change notification settings - Fork 98
/
Copy pathlyapunov_symmetric_eig.m
41 lines (34 loc) · 1.26 KB
/
lyapunov_symmetric_eig.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
function X = lyapunov_symmetric_eig(V, lambda, C, tol)
% 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
% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 31, 2018.
% Contributors:
% Change log:
% Sep. 24, 2023 (NB):
% Edited out bsxfun().
% AX + XA = C is equivalent to DY + YD = M with
% Y = V'XV, M = V'CV and D = diag(lambda).
M = V'*C*V;
% W(i, j) = lambda(i) + lambda(j)
W = lambda + lambda';
% Normally, the solution Y is simply this:
Y = M ./ W;
% But this may involve divisions by (almost) 0 in certain places.
% Thus, we go for a pseudo-inverse.
absW = abs(W);
if ~exist('tol', 'var') || isempty(tol)
tol = numel(C)*eps(max(absW(:))); % similar to pinv tolerance
end
Y(absW <= tol) = 0;
% Undo the change of variable
X = V*Y*V';
end