## 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://software.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 |

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 `I`_{k} 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 `X``B` - `A``X` = `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 `C`^{T}`C` + `S`^{T}`S` = `I`_{N-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 `α``A``X` + `β``X``B` = `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 `M``N` linear equations for `M``N` unknown components of the vector `x` and right hand side vector `f`:

`x` = (`x`_{11}, `x`_{21}, ..., `x`_{M1}, `x`_{12}, `x`_{22}, ..., `x`_{M2}, ..., `x`_{1N}, `x`_{2N}, ..., `x`_{MN})^{T}

`f` = (`f`_{11}, `f`_{21}, ..., `f`_{M1}, `f`_{12}, `f`_{22}, ..., `f`_{M2}, ..., `f`_{1N}, `f`_{2N}, ..., `f`_{MN})^{T}

Matrix of order `M``N` 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 `I`_{M} represents the unit matrix of order `M`. Thus the coefficient matrix is

This matrix is sparse, with `M` + `N` - 1 nonzero elements in each `M``N`-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