Home > manopt > tools > lyapunov_symmetric.m

# lyapunov_symmetric

## PURPOSE Solves AX + XA = C when A = A', as a pseudo-inverse.

## SYNOPSIS function X = lyapunov_symmetric(A, C, tol)

## DESCRIPTION ``` 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.

## CROSS-REFERENCE INFORMATION This function calls:
This function is called by:

## SOURCE CODE ```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 %
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```

Generated on Tue 19-May-2020 18:46:12 by m2html © 2005