Solves AX + XA = C when A = A', as a pseudo-inverse. function X = lyapunov_symmetric(A, C) function X = lyapunov_symmetric(A, C, tol) Matrices A, C and X have size nxn. When the solution exists and is unique, this is equivalent to sylvester(A, A', C) or lyap(A, -C). This works for both real and complex inputs. If C is a 3-D array of size nxnxk, then X has size nxnxk as well, and each slice X(:, :, i) corresponds to the solution for the system with right-hand side C(:, :, i). This is more efficient then calling the function multiple times for each slice of C. If the solution is not unique, the smallest-norm solution is returned. If a solution does not exist, a minimum-residue solution is returned. tol is a tolerance used to determine which eigenvalues are numerically zero. It can be specified manually; otherwise, a default value is chosen. Overall, if A is nxn, the output is very close to: X = reshape(pinv(kron(A, eye(n)) + kron(eye(n), A))*C(:), [n, n]), but it is computed far more efficiently. The most expensive step is an eigendecomposition of A, whose complexity is essentially O(n^3) flops. If A is not symmetric, only its symmetric part is used: (A+A')/2. If C is (skew-)symmetric, then X is (skew-)symmetric (up to round-off), and similarly in the complex case. To solve one system at a time, while reusing the eigendecomposition of A, call lyapunov_symmetric_eig. See also: lyapunov_symmetric_eig sylvester lyap sylvester_nochecks
0001 function X = lyapunov_symmetric(A, C, tol) 0002 % Solves AX + XA = C when A = A', as a pseudo-inverse. 0003 % 0004 % function X = lyapunov_symmetric(A, C) 0005 % function X = lyapunov_symmetric(A, C, tol) 0006 % 0007 % Matrices A, C and X have size nxn. When the solution exists and is 0008 % unique, this is equivalent to sylvester(A, A', C) or lyap(A, -C). 0009 % This works for both real and complex inputs. 0010 % 0011 % If C is a 3-D array of size nxnxk, then X has size nxnxk as well, and 0012 % each slice X(:, :, i) corresponds to the solution for the system with 0013 % right-hand side C(:, :, i). This is more efficient then calling the 0014 % function multiple times for each slice of C. 0015 % 0016 % If the solution is not unique, the smallest-norm solution is returned. 0017 % 0018 % If a solution does not exist, a minimum-residue solution is returned. 0019 % 0020 % tol is a tolerance used to determine which eigenvalues are numerically 0021 % zero. It can be specified manually; otherwise, a default value is chosen. 0022 % 0023 % Overall, if A is nxn, the output is very close to: 0024 % X = reshape(pinv(kron(A, eye(n)) + kron(eye(n), A))*C(:), [n, n]), 0025 % but it is computed far more efficiently. The most expensive step is an 0026 % eigendecomposition of A, whose complexity is essentially O(n^3) flops. 0027 % 0028 % If A is not symmetric, only its symmetric part is used: (A+A')/2. 0029 % 0030 % If C is (skew-)symmetric, then X is (skew-)symmetric (up to round-off), 0031 % and similarly in the complex case. 0032 % 0033 % To solve one system at a time, while reusing the eigendecomposition of A, 0034 % call lyapunov_symmetric_eig. 0035 % 0036 % See also: lyapunov_symmetric_eig sylvester lyap sylvester_nochecks 0037 0038 % This file is part of Manopt: www.manopt.org. 0039 % Original author: Nicolas Boumal, April 17, 2018. 0040 % Contributors: 0041 % Change log: 0042 % Aug. 31, 2018 (NB): 0043 % Now works with C having multiple slices (nxnxk), and added some 0044 % safeguards. 0045 0046 n = size(A, 1); 0047 assert(ismatrix(A) && size(A, 2) == n, 'A must be square.'); 0048 assert(size(C, 1) == n && size(C, 2) == n, ... 0049 'Each slice of C must have the same size as A.'); 0050 0051 if ~exist('tol', 'var') 0052 tol = []; 0053 end 0054 0055 % Make sure A is numerically Hermitian (or symmetric). 0056 % The cost of this safety step is negligible compared to eig. 0057 A = (A+A')/2; 0058 0059 % V is unitary or orthogonal and lambda is real. 0060 [V, lambda] = eig(A, 'vector'); 0061 0062 % Solve for each slice separately. 0063 X = zeros(size(C)); 0064 for k = 1 : size(C, 3) 0065 X(:, :, k) = lyapunov_symmetric_eig(V, lambda, C(:, :, k), tol); 0066 end 0067 0068 end 0069