Computing principal angles between invariant subspaces of block triangular matrices
Goal
Get information about the relative position of two invariant subspaces of a block triangular matrix.
Solution
Assuming the subspaces are represented as spans of some vectors, the relative position of the subspaces can be obtained via calculating the set of principal angles between the subspaces (see Computing principal angles between two subspaces). Additionally, for invariant subspaces of block triangular matrices the Sylvester matrix equation must be solved. The solver used depends on the matrix characteristics:
If both diagonal blocks of the triangular matrix are upper triangular, use the LAPACK ?trsyl routine.
If both diagonal blocks of the triangular matrix are not large and not upper triangular, use LAPACK linear solvers.
If both diagonal blocks of the triangular matrix are large, upper triangular, and sparse, use the oneMKL PARDISO solver.
Source code: see these files in the samples archive available at https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip:
ANGLES/uep_subspace1/main.f
ANGLES/uep_subspace2/main.f
ANGLES/uep_subspace3/main.f
Solving Sylvester matrix equation using LAPACK ?trsyl
CALL DTRSYL('N', 'N', -1, K, N-K, AA, N, AA(K+1,K+1), N,
& AA(1,K+1), N, ALPHA, INFO)
IF(INFO.EQ.0) THEN
PRINT *,"DTRSYL completed, SCALE=",ALPHA
ELSE IF(INFO.EQ.1) THEN
PRINT *,"DTRSYL solved perturbed equations"
ELSE
PRINT *,"DTRSYL failed. INFO=",INFO
STOP
END IF
Solving Sylvester matrix equation using LAPACK linear solvers
REAL*8 AA(N,N), FF(K*(N-K)), AAA(K*(N-K),K*(N-K))
…
! Forming dense coefficient matrix for Sylvester equation
CALL SYLMAT(K, AA, N, N-K, AA(K+1,K+1), N, -1D0, 1D0, AAA, NK,
& INFO)
! Processing INFO returned by SYLMAT
…
! Forming the right hand side for the system of linear equations that
! correspond to Sylvester equation.
DO I = 1, K
DO J = 1, N-K
FF((J-1)*K+I) = AA(I,J+K)
END DO
END DO
! Solve the system of linear equations
CALL DGESV(NK, 1, AAA, NK, IPIV, FF, NK, INFO )
! Processing INFO returned by DGESV
Solving Sylvester matrix equation using oneMKL PARDISO
REAL*8 AA(N,N), FF(K*(N-K)), VAL(K*(N-K)*(N-1))
INTEGER ROWINDEX(K*(N-K)+1), COLS(K*(N-K)*(N-1))
…
! Forming sparse coefficient matrix for Sylvester equation
CALL FSYLVOP(K, AA, N, N-K, AA(K+1,K+1), N, -1D0, 1D0, COLS,
& ROWINDEX, VAL, INFO)
! Processing INFO returned by FSYLVOP
…
! Form the right hand side of the Sylvester equation
DO I=1,K
DO J=1,N-K
FF((J-1)*K+I) = AA(I,J+K)
END DO
END DO
CALL PARDISOINIT (PT, 1, IPARM)
CALL PARDISO (PT, 1, 1, 11, 13, NK, VAL, ROWINDEX,
& COLS, PERM, 1, IPARM, 1, FF, X, IERR)
! Processing IERR returned by PARDISO
…
Discussion
Task |
Routine |
Description |
|---|---|---|
Solve Sylvester matrix equation for matrix with upper triangular diagonal blocks |
dtrsyl |
Solve Sylvester equation for real quasi-triangular or complex triangular matrices. |
Solve Sylvester matrix equation for matrix which is small and not upper triangular |
dgesv |
Computes the solution to the system of linear equations with a square matrix A and multiple right-hand sides. |
Solve Sylvester matrix equation for matrix which is not small and not upper triangular |
pardiso |
Calculates the solution of a set of sparse linear equations with single or multiple right-hand sides. |
In order to determine the principal angles between invariant subspaces of the matrix, first let an N-by-N matrix be represented in a block triangular form:
Here diagonal blocks A and B are square matrices of order k and N-k, respectively. If Ik denotes the unit matrix of order k, the equality
means the span of first k vectors of the standard basis is invariant with respect to transformations on matrix A.
Another invariant subspace can be found as a span of columns of the compound matrix 
Here X is some rectangular k-by-(N - k) matrix which should be found. Compute the product
If X is a solution of the Sylvester equation XB - AX = F the result in the last equation is
This demonstrates invariance of the subspace spanned by columns of 
QR factorization can be used to orthogonalize the basis in the second invariant subspace:
Here C is a k-by-(N-k) matrix and S is an (N-k)-by-(N-k) matrix. C and S satisfy the equation CTC + STS = IN-k, where R is an upper triangular square matrix of order N-k. Compute principal angles between these two invariant subspaces using the SVD of C:
Diagonal elements of Σ are the cosines of the principal angles.
Matrix of Sylvester equations
Consider the Sylvester equation αAX + βXB = F.
Here square matrices A and B have orders M and N, respectively, and α and β are scalars. F is a given M-by-N matrix:
X is the M-by-N matrix to be found:
This matrix equation can be considered as a system 
of MN linear equations for MN unknown components of the vector x and right hand side vector f:
x = (x11, x21, ..., xM1, x12, x22, ..., xM2, ..., x1N, x2N, ..., xMN)T
f = (f11, f21, ..., fM1, f12, f22, ..., fM2, ..., f1N, f2N, ..., fMN)T
Matrix of order MN can be represented as a sum of two matrices. One corresponds to multiplication of matrix X from the left by matrix A which can be represented in block form with blocks of size M by M. The matrix forms a block diagonal matrix with N blocks on the diagonal:
The other matrix in the sum corresponds to multiplication of matrix X from the right by matrix B. Using the same block form representation yields
Here IM represents the unit matrix of order M. Thus the coefficient matrix is
This matrix is sparse, with M + N - 1 nonzero elements in each MN-element row. Therefore the oneMKL PARDISO sparse solver can be used effectively (see the code ANGLES/source/fsylvop.f, which forms the coefficient matrix in CSR format for oneMKL PARDISO). However, for comparatively small M and N the oneMKL LAPACK linear solvers are more efficient (see the code ANGLES/source/sylmat.f, which forms the coefficient matrix as a dense matrix for use with dgesv).
Parent topic: Intel® oneAPI Math Kernel Library Cookbook