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 See also: stiefelfactory rotationsfactory sylvester sylvester_nochecks
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 % 0022 % See also: stiefelfactory rotationsfactory sylvester sylvester_nochecks 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