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

Routines Used

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
Level Two Title