lyapunov_symmetric.m

function X = lyapunov_symmetric(A, C, tol)
% 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

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 17, 2018.
% Contributors: 
% Change log: 
%   Aug. 31, 2018 (NB):
%       Now works with C having multiple slices (nxnxk), and added some
%       safeguards.

    n = size(A, 1);
    assert(ismatrix(A) && size(A, 2) == n, 'A must be square.');
    assert(size(C, 1) == n && size(C, 2) == n, ...
           'Each slice of C must have the same size as A.');
       
   if ~exist('tol', 'var')
       tol = [];
   end

    % Make sure A is numerically Hermitian (or symmetric).
    % The cost of this safety step is negligible compared to eig.
    A = (A+A')/2;

    % V is unitary or orthogonal and lambda is real.
    [V, lambda] = eig(A, 'vector');
    
    % Solve for each slice separately.
    X = zeros(size(C));
    for k = 1 : size(C, 3)
        X(:, :, k) = lyapunov_symmetric_eig(V, lambda, C(:, :, k), tol);
    end

end