Home > manopt > manifolds > stiefel > solve_for_triu.m

# solve_for_triu

## PURPOSE

Solve the linear matrix equation AX + X'A' = H for X upper triangular.

## SYNOPSIS

function X = solve_for_triu(A, H)

## DESCRIPTION

``` Solve the linear matrix equation AX + X'A' = H for X upper triangular.

function X = solve_for_triu(A, H)

Given A of size p-by-p and H (symmetric) of size p-by-p, solves the
linear matrix equation AX + X'A' = H for X upper triangular.

The total computational cost is O(p^4).

If the same equation is to be solved for X symmetric instead, call
Matlab's built-in sylvester function.

This is a support function to compute the inverse of QR-based
retractions.

This algorithm appears as Algorithm 1 in:
Empirical Arithmetic Averaging over the Compact Stiefel Manifold,
Tetsuya Kaneko, Simone Fiori, Toshihisa Tanaka,
IEEE Transactions on Signal Processing, 2013

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:
• rotationsfactory Returns a manifold structure to optimize over rotation matrices.
• stiefelfactory Returns a manifold structure to optimize over orthonormal matrices.

## SOURCE CODE

```0001 function X = solve_for_triu(A, H)
0002 % Solve the linear matrix equation AX + X'A' = H for X upper triangular.
0003 %
0004 % function X = solve_for_triu(A, H)
0005 %
0006 % Given A of size p-by-p and H (symmetric) of size p-by-p, solves the
0007 % linear matrix equation AX + X'A' = H for X upper triangular.
0008 %
0009 % The total computational cost is O(p^4).
0010 %
0011 % If the same equation is to be solved for X symmetric instead, call
0012 % Matlab's built-in sylvester function.
0013 %
0014 % This is a support function to compute the inverse of QR-based
0015 % retractions.
0016 %
0017 % This algorithm appears as Algorithm 1 in:
0018 %  Empirical Arithmetic Averaging over the Compact Stiefel Manifold,
0019 %  Tetsuya Kaneko, Simone Fiori, Toshihisa Tanaka,
0020 %  IEEE Transactions on Signal Processing, 2013
0021 %
0023
0024 % This file is part of Manopt: www.manopt.org.
0025 % Original author: Nicolas Boumal, July 18, 2018.
0026 % Contributors:
0027 % Change log:
0028 %
0029 %   Aug. 3, 2018 (NB):
0030 %       Initial array of zeros now copies the class of A, so that if A is a
0031 %       regular matrix of doubles it doesn't change anything, but if A is
0032 %       on GPU, then this function also works on GPU.
0033
0034 % One tentative idea to reduce the cost to O(p^3) involves taking an LU
0035 % factorization of A. Then, we obtain a permutation matrix P and
0036 % triangular matrices L (lower) and U (upper) such that PA = LU.
0037 % Since inv(P) = P', the matrix equation becomes:
0038 %
0039 %  P' L U X + X' U' L' P = H
0040 %
0041 % Notice that U*X is still upper triangular, so that we can solve for
0042 % U*X first, and obtain X later by solving an upper triangular system.
0043 % After this change of variables, the system involves P'L instead of A.
0044 % If the permutation happens to be identity, then clearly principal
0045 % submatrices of P'L = L are lower triangular themselves, and as a
0046 % result the linear systems that we need to solve below only cost
0047 % O(pp^2) instead of O(pp^3). Summing for pp = 1 .. p gives O(p^3)
0048 % instead of O(p^4). In general though, P is not the identity
0049 % permutation. In particular, if P' reverses the order of the rows of L,
0050 % so that the first half of the principal submatrices of P'L are full,
0051 % then we revert back to O(p^4) complexity. Interestingly, for X, Y close
0052 % by, the matrix A = X'*Y that appears in computing the inverse retraction
0053 % is close to identity, so that its LU factorization naturally leads to P
0054 % identity; thus, in such scenario we could reduce the cost to O(p^3)
0055 % without loss of stability due to the LU change of variable.
0056
0057     p = size(A, 1);
0058     X = zeros(p, p, class(A));
0059     for pp = 1 : p
0060         b = H(1:pp, pp);
0061         b(end) = b(end)/2;
0062         b(1:end-1) = b(1:end-1) - X(1:pp-1, 1:pp-1)'*A(pp, 1:pp-1)';
0063         X(1:pp, pp) = A(1:pp, 1:pp) \ b;
0064     end
0065
0066 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005