Intel® oneAPI Math Kernel Library Cookbook

Download PDF
ID 758503
Date 3/31/2023
Public

A newer version of this document is available. Customers should click here to go to the newest version.

Document Table of Contents
Document Table of Contents x
Intel® oneAPI Math Kernel Library Cookbook
Intel® oneAPI Math Kernel Library Cookbook x
Notices and Disclaimers Getting Technical Support What's New Notational Conventions Related Information Finding an approximate solution to a stationary nonlinear heat equation Factoring general block tridiagonal matrices Solving a system of linear equations with an LU-factored block tridiagonal coefficient matrix Factoring block tridiagonal symmetric positive definite matrices Solving a system of linear equations with a block tridiagonal symmetric positive definite coefficient matrix Computing principal angles between two subspaces Computing principal angles between invariant subspaces of block triangular matrices Evaluating a Fourier integral Using Fast Fourier Transforms for computer tomography image reconstruction Noise filtering in financial market data streams Using the Monte Carlo method for simulating European options pricing Using the Black-Scholes formula for European options pricing Multiple simple random sampling without replacement Using a histospline technique to scale images Speeding up Python* scientific computations Bibliography
Intel® oneAPI Math Kernel Library Cookbook

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:

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

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

    2. Using the threshold, determine the dimensions of the subspaces.

    3. Form the orthonormal bases.

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

  3. 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

Routines Used

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 UP = QR. 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 UP = QR means that all columns of U are linear combinations of the first l columns of Q. Due to pivoting, the diagonal elements rj,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 |rl,l| > threshold > |rl+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 u1 and w1 are called principal vectors. The other principal angles and vectors are defined recursively:

The principal vectors from the same subspace are pairwise orthogonal:

(ui, uj) = (wi, wj) = δ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 UTW:

UTW = PΣQT, PTP = Ik, QQT = Il

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


The respective pairs of principal vectors are Upi and Wqi where pi and qi are the i-th columns of P and Q.

Parent topic: Intel® oneAPI Math Kernel Library Cookbook
Level Two Title