## Computing principal angles between two subspaces

### Goal

Get information about the relative position of two subspaces in an inner product space.

### Solution

Assuming the subspaces are represented as spans of some vectors, the relative position of the subspaces can be obtained by calculating the set of Principal Angles between the subspaces. To calculate the angles:

Build orthonormal bases in each subspace and determine the dimensions of the subspaces.

Call an appropriate subroutine to perform a QR factorization with pivoting of matrices, the columns of which span the subspaces.

Using the threshold, determine the dimensions of the subspaces.

Form the orthonormal bases.

Form a matrix of inner products of the basis vectors from the one subspace and the basis vectors of another subspace.

Compute the Singular Value Decomposition of the matrix.

Source code: see the ANGLES/definition/main.f file in the samples archive available at https://software.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip.

### Building orthonormal bases and determining subspace dimensions

```
…
REAL*8 Y(LDY,K),TAU(N),WORK(3*N)
…
!
! Apply QR factorization with pivoting
CALL DGEQPF(N, K, Y, LDY, JPVT, TAU, WORK, INFO)
! Process info returned by DGEQPF
…
! Compute the rank
K1=0
DO WHILE((K1 .LT. K) .AND. (ABS(Y(K1+1,K1+1)) .GT. THRESH))
K1 = K1 + 1
END DO
!
! Form K1 orthonormal vectors via call DORGQR
CALL DORGQR(N, K1, K1, Y, LDY, TAU, WORK, LWORK, INFO)
! Process info returned by DORGQR
…
```

### Forming matrix of inner products and computing SVD

```
REAL*8 U(N,KU),V(N,KV),W(KU,KV),VECL(KU,KMIN)
REAL*8 VECRT(KMIN,KV),S(KMIN),WORK(5*KU)
…
! Form W=U^t*V
CALL DGEMM(‘T’, ‘N’, KU, KV, N, 1D0, U, N, V, N, 0D0, W, KU1)
…
! SVD of W=U^t*V
CALL DGESVD(‘S’, ‘S’, KU, KV, W, KU, S, VECL, KU, VECRT, KMIN, WORK, LWORK, INFO)
! Process info returned by DGESVD
…
```

### Discussion

Task |
Routine |
Description |
---|---|---|

QR factorization with pivoting of matrices |
dgeqpf |
Compute the QR factorization of a general m-by-n matrix with pivoting |

Form orthonormal bases |
dorgqr |
Generates the real orthogonal matrix Q of the QR factorization formed by ?geqpf or ?geqp3 |

Form a matrix of inner products of the basis vectors from the one subspace and the basis vectors of another subspace. |
dgemm |
Compute a matrix-matrix product with general matrices |

Compute the Singular Value Decomposition of the matrix |
dgesvd |
Compute the singular value decomposition of a general rectangular matrix |

The first step is to build orthonormal bases in each subspace and determine the dimensions of the subspaces.

Let `U` be an `N`-by-`k` matrix (`N`≥`k`) with columns representing vectors in some inner product linear space. To construct an orthonormal basis in this space you can use QR factorization of the matrix `U`, which with pivoting can be represented as `U``P` = `Q``R`. If the dimension of the space is `l` (`l`≤`k`), in the absence of rounding errors occur this yields an orthogonal (unitary for complex-valued matrices) `N`-by-`N` matrix `Q` and upper triangular `N`-by-`k` matrix `R`:

The equation `U``P` = `Q``R` means that all columns of `U` are linear combinations of the first `l` columns of `Q`. Due to pivoting, the diagonal elements `r`_{j,j} of `R` are ordered in non-increasing order of absolute values. In fact, pivoting provides even stronger inequalities:

for `j`≤`m`≤`k`.

In actual computations with rounding errors, the elements of the lower right (`k` - `l `)-by-(`k` - `l `) triangle of `R` are small but non-zero, so a threshold is used to determine the rank |`r`_{l,l}| > `threshold` > |`r`_{l+1,l+1}|.

Now you can determine the set of angles between subspaces.

Let and be two subspaces in the same `N`-dimensional Euclidean space with dim()=`k`, dim()=`l`, and `k`≤`l`. To find out the relative position of these subspaces you can use principal angles `θ`_{1}≥`θ`_{2}≥...≥`θ`_{k}≥ 0, which are defined as follows.

The first angle is defined as:

The vectors `u`_{1} and `w`_{1} are called principal vectors. The other principal angles and vectors are defined recursively:

The principal vectors from the same subspace are pairwise orthogonal:

(`u`_{i}, `u`_{j}) = (`w`_{i}, `w`_{j}) = `δ`_{ij}

To compute the principal angles you can use Singular Value Decomposition of matrices. Let `U` and `W` be matrices of sizes `N`-by-`k` and `N`-by-`l` respectively, with columns being orthonormal bases in and respectively. Compute the SVD of the `k`-by-`l` matrix `U`^{T}`W`:

`U`^{T}`W` = `P``Σ``Q`^{T}, `P`^{T}`P` = `I`_{k}, `Q``Q`^{T} = `I`_{l}

It can be proven that the diagonal elements of `Σ` are the cosines of the principal angles:

The respective pairs of principal vectors are `U``p`^{i} and `W``q`^{i} where `p`^{i} and `q`^{i} are the `i`-th columns of `P` and `Q`.

**Parent topic:**Intel® oneAPI Math Kernel Library Cookbook