Next Article in Journal
(ρ,η,μ)-Interpolative Kannan Contractions I
Next Article in Special Issue
Using Free Mathematical Software in Engineering Classes
Previous Article in Journal
Qualitative and Quantitative Analyses of COVID-19 Dynamics
Previous Article in Special Issue
Gauss–Newton–Secant Method for Solving Nonlinear Least Squares Problems under Generalized Lipschitz Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Algorithms for Computing an Arbitrary Singular Value of a Tensor Sum †

1
National Institute of Technology, Kagawa College, Takuma Campus, Mitoyo 769-1192, Japan
2
Department of Applied Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan
*
Author to whom correspondence should be addressed.
This paper is an extended version of our paper published in International Conference on Mathematics and Its Applications in Science and Engineering (ICMASE2020), Ankara Haci Bayram Veli University, Turkey (Online), 9–10 July 2020.
Axioms 2021, 10(3), 211; https://doi.org/10.3390/axioms10030211
Submission received: 19 July 2021 / Revised: 22 August 2021 / Accepted: 26 August 2021 / Published: 31 August 2021
(This article belongs to the Special Issue Numerical Analysis and Computational Mathematics)

Abstract

:
We consider computing an arbitrary singular value of a tensor sum: T : = I n I m A + I n B I + C I m I R m n × m n , where A R × , B R m × m , C R n × n . We focus on the shift-and-invert Lanczos method, which solves a shift-and-invert eigenvalue problem of ( T T T σ ˜ 2 I m n ) 1 , where σ ˜ is set to a scalar value close to the desired singular value. The desired singular value is computed by the maximum eigenvalue of the eigenvalue problem. This shift-and-invert Lanczos method needs to solve large-scale linear systems with the coefficient matrix T T T σ ˜ 2 I m n . The preconditioned conjugate gradient (PCG) method is applied since the direct methods cannot be applied due to the nonzero structure of the coefficient matrix. However, it is difficult in terms of memory requirements to simply implement the shift-and-invert Lanczos and the PCG methods since the size of T grows rapidly by the sizes of A, B, and C. In this paper, we present the following two techniques: (1) efficient implementations of the shift-and-invert Lanczos method for the eigenvalue problem of T T T and the PCG method for T T T σ ˜ 2 I m n using three-dimensional arrays (third-order tensors) and the n-mode products, and (2) preconditioning matrices of the PCG method based on the eigenvalue and the Schur decomposition of T. Finally, we show the effectiveness of the proposed methods through numerical experiments.

1. Introduction

We consider computing an arbitrary singular value of a tensor sum:
T : = I n I m A + I n B I + C I m I R m n × m n ,
where A R × , B R m × m , C R n × n , I n is the n × n identity matrix, and the symbol “⊗” denotes the Kronecker product. The tensor sum T arises from a finite difference discretization of three-dimensional constant coefficient partial differential equations (PDE) defined as follows:
[ a · ) + b · + c u ( x , y , z ) = g ( x , y , z ) in Ω , u ( x , y , z ) = 0 on Ω ,
where Ω = ( 0 , 1 ) × ( 0 , 1 ) × ( 0 , 1 ) , a , b R 3 , c R , and the symbol “*” denotes elementwise products. If a = ( 1 , 1 , 1 ) , then a · ( ) = Δ . Matrix T tends to be too large even if A , B and C are not. Hence it is difficult to compute singular values of T with regard to the memory requirement.
Previous studies [1,2] provided methods to compute the maximum and minimum singular values of T. By the previous studies, one can compute only the maximum and minimum singular values of T without shift. On the other hand, one can compute arbitrary singular values of T with the shift by this work. The previous studies are based on the Lanczos bidiagonalization method (see, e.g., [3]), which computes the maximum and minimum singular values of a matrix. For insights on Lanczos bidiagonalization method, see, e.g., [4,5,6]. The Lanczos bidiagonalization method for T was implemented using tensors and their operations to reduce the memory requirement.
The Lanczos method with the shift-and-invert technique, see, e.g., [3], is widely known for computing an arbitrary eigenvalue λ of a symmetric matrix M R n × n . This method solves the shift-and-invert eigenvalue problem: ( M σ ˜ I n ) 1 x = ( λ σ ˜ ) 1 x , where x is the eigenvector of M corresponding to λ , and σ ˜ is a shift point which is set to the nearby λ ( σ ˜ λ ). Since the eigenvalue problem has the eigenvalue ( λ σ ˜ ) 1 as the maximum eigenvalue, the method is effective for computing the desired eigenvalue λ near σ ˜ . For successful work using the shift-and-invert technique, see, e.g., [7,8,9,10,11,12,13].
Therefore, we obtain a computing method for an arbitrary singular value of T based on the shift-and-invert Lanczos method. The method solves the following shift-and-invert eigenvalue problem: ( T T T σ ˜ 2 I m n ) 1 x = ( σ 2 σ ˜ 2 ) 1 x , where σ is the desired singular value of T, x is the corresponding right-singular vector, and σ ˜ is close to σ ( σ ˜ σ ). This shift-and-invert Lanczos method requires the solution of large-scale linear systems with the coefficient matrix T T T σ ˜ 2 I m n . Here, T T T σ ˜ 2 I m n can be a dense matrix whose number of elements is O ( n 6 ) even if T is a sparse matrix whose number of elements is O ( n 4 ) when A , B , C R n × n are dense.
Since it is difficult regarding the memory requirement to apply the direct method, e.g., the Cholesky decomposition, which needs generating matrix T T T σ ˜ 2 I m n , the preconditioned conjugate gradient (PCG) method, see, e.g., [14], is applied, even though it is difficult in terms of memory requirements to simply implement this shift-and-invert Lanczos method and the PCG method since the size of T grows rapidly by the sizes of A, B, and C.
We propose the following two techniques in this paper: (1) Efficient implementations of the shift-and-invert Lanczos method for the eigenvalue problem of T T T and the PCG method for T T T σ ˜ 2 I m n using three-dimensional arrays (third-order tensors) and the n-mode products, see, e.g., [15]. (2) Preconditioning matrices based on the eigenvalue decomposition and the Schur decomposition of T for faster convergence of the PCG method. Finally, we show the effectiveness of the proposed method through numerical experiments.

2. Preliminaries of Tensor Operations

A tensor means a multidimensional array. Particularly, the third-order tensor X R I × J × K plays an important role. In the rest of this section, the definitions of some tensor operations are shown. For more details, see, e.g., [15].
Firstly, a summation, a subtraction, an inner product, and a norm for X , Y R I × J × K are defined as follows:
( X ± Y ) i j k : = X i j k ± Y i j k , ( X , Y ) : = i = 1 I j = 1 J k = 1 K X i j k Y i j k , X = ( X , X ) ,
where X i j k denotes the ( i , j , k ) element of X . Secondly, the n-mode product of a tensor X R I 1 × I 2 × × I N and a matrix M R J × I n is defined as
( X × n M ) i 1 i n 1 j i n + 1 i N = i n = 1 I n X i 1 i 2 i N M p i n ,
where n { 1 , 2 , N } , i k { 1 , 2 , , I k } for k = 1 , 2 , N , and j { 1 , 2 , , J } . Finally, vec and vec 1 operators are the following maps between a vector space R I J K and a tensor space R I × J × K : vec : R I × J × K R I J K and vec 1 : R I J K R I × J × K .   vec operator can vectorize a tensor by combining all column vectors of the tensor into one long vector. Conversely, vec 1 operator can reshape a tensor from one long vector.

3. Shift-and-Invert Lanczos Method for an Arbitrary Singular Value over Tensor Space

This section gives an algorithm for computing an arbitrary singular value of the tensor sum T. Let σ and x be a desired singular value of T and the corresponding right singular vectors, respectively. Then, the eigenvalue problem of T is written by T T T x = σ 2 x . Here, introducing a shift σ ˜ σ , the shift-and-invert eigenvalue problem is
( T T T σ ˜ 2 I m n ) 1 x = 1 σ 2 σ ˜ 2 x .
The shift-and-invert Lanczos method (see, e.g., [3]) computes the nearest singular value σ based on Equation (3). Reconstructing this method over the × m × n tensor space, we obtain Algorithm 1 whose memory requirement is of O ( n 3 ) when n = m = .
Algorithm 1: Shift-and-invert Lanczos method for an arbitrary singular value over tensor space
1:  Choose an initial tensor Q 0 R × m × n ;
2:  V : = Q 0 , β 0 : = | | V | | ;
3:  for k = 1 , 2 , , until convergence do
4:   Q k : = V / β k 1 ;
5:   V : = vec 1 T T T σ ˜ 2 I m n 1 vec ( Q k ) ;
     (Computed by Algorithms 3 or 4 in Section 4)
6:   V : = V β k 1 Q k 1 ;
7:   α k : = ( Q k , V ) ;
8:   V : = V α k Q k ;
9:   β k : = | | V | | ;
10:  end for
11:  Approximate singular value σ = σ ˜ 2 + 1 λ ˜ ( k ) , where λ ˜ ( k ) is the maximum eigenvalue of T ˜ k .
At step k, we have the following T ˜ k by Algorithm 1:
T ˜ k : = α 1 β 1 β 1 α 2 β 2 β k 2 α k 1 β k 1 β k 1 α k R k × k .
To implement Algorithm 1, we need to iteratively solve the linear system
V : = vec 1 T T T σ ˜ I m n 1 vec ( Q k ) ,
whose coefficient matrix is m n × m n , that is, the memory requirement is O ( n 6 ) when n = m = . Here, the convergence rate of the shift and invert Lanczos method depends on the ratio of gaps between the maximum, the second maximum, and the minimum singular values σ 1 , σ 2 , σ m of ( T T T σ ˜ 2 I m n ) 1 as follows: ( σ 1 2 σ 2 2 ) / ( σ 1 2 σ m 2 ) .
In the next section, we consider solving the linear systems with memory requirement of O ( n 3 ) when n = m = .

4. Preconditioned Conjugate Gradient (PCG) Method over Tensor Space

This section provides an efficient solver of Equation (4) using tensors. This linear system is rewritten by v = T T T σ ˜ 2 I m n 1 q k , where v : = vec ( V ) and q k : = vec ( Q k ) . Then we solve T T T σ ˜ 2 I m n v = q k , where v and q k are unknown and known vectors. Since the coefficient matrix is symmetric positive definite, we can use the preconditioned conjugate gradient method (PCG method, see, e.g., [14]), which is one of the widely used solvers. However, it is difficult to simply apply the method due to the complex nonzero structure of the coefficient matrix T T T σ ˜ 2 I m n . For applying the PCG method, we consider transforming the linear system T T T σ ˜ 2 I m n v = q k by the eigendecomposition and the complex Schur decomposition as shown in the next subsections.

4.1. PCG Method by the Eigendecomposition

Firstly, T is decomposed into T : = X D X 1 , where X and D are a matrix whose column vectors are eigenvectors and a diagonal matrix with eigenvalues, respectively. Then, it follows that
T T T σ ˜ 2 I m n v = q k ( X D X 1 ) H ( X D X 1 ) σ ˜ 2 I m n v = q k D ¯ X H X D σ ˜ 2 X H X X 1 v = X H q k ,
where D ¯ is the complex conjugate of D. We rewrite the above linear system into A ˜ y ˜ = b ˜ , where A ˜ : = D ¯ X H X D σ ˜ 2 X H X , y ˜ : = X 1 v , and b ˜ : = X H q k . Here, X is easily computed by small matrices X A , X B , and X C whose column vectors are eigenvectors of A, B, and C as follows: X = X C X B X A . Moreover, eigenvalues of T in D are obtained by summations of each eigenvalue of A, B, and C.
The PCG method for solving A ˜ y ˜ = b ˜ is shown in Algorithm 2. Since this algorithm computes y ˜ , we need to compute v = X y ˜ . Section 4.1.1 proposes a preconditioning matrix and Section 4.1.2 provides efficient computations using tensors.
Algorithm 2: PCG method over vector space for A ˜ y ˜ = b ˜
1:  Choose an initial vector x 0 R m n and p 0 = 0 R m n , and an initial scalar β 0 = 0 ;
2:  r 0 = b ˜ A ˜ x 0 ;
3:  z 0 = M 1 r 0 ;
4:  for k = 1 , 2 , , until convergence do
5:   p k = z k 1 + β k 1 p k 1 ;
6:   p ^ k = A ˜ p k ;
7:   α k = z k 1 , r k 1 / p k 1 , p ^ k ;
8:   x k = x k 1 + α k p k ;
9:   r k = r k 1 α k p ^ k ;
10:   z k = M 1 r k ;
11:   β k = z k , r k / z k 1 , r k 1 ;
12:  end for
13:  Obtain an approximate solution y ˜ x k ;

4.1.1. Preconditioning Matrix

Algorithm 2 solves
M 1 A ˜ M H M H y ˜ = M 1 b ˜ ,
where M R m n × m n is a preconditioning matrix. M must satisfy the following two conditions: (1) a condition number of M 1 A ˜ is close to 1; (2) the matrix-vector multiplication of M 1 is easily computed.
Therefore, we propose a preconditioning matrix based on the eigendecomposition of T
M : = D ¯ D σ ˜ 2 I m n .
Since M is the diagonal matrix, the second condition of the preconditioning matrix is satisfied. Moreover, if T is symmetric, X is the unitary matrix, that is, X H X = I m n . In the case of the symmetric matrix T, we obtain M = A ˜ . Namely, the proposed matrix satisfies the first conditions when T is symmetric. So, even if T is not exactly symmetric, if T is almost symmetric, then we can expect the preconditioning matrix M to be effective.

4.1.2. Efficient Implementation of Algorithm 2 by the Eigendecomposition

Similarly to obtaining Algorithm 1, to improve an implementation of Algorithm 2, we reconstruct m n dimensional vectors into × m × n tensors via vec 1 operator as follows: X k : = vec 1 ( x k ) , R k : = vec 1 ( r k ) , P k : = vec 1 ( p k ) , Z k : = vec 1 ( z k ) , and P ^ k : = vec 1 ( p ^ k ) . Most computations of vectors are simply transformed into computations of tensors because of the linearity of vec 1 operator.
In the rest of this section, we show the computations of vec 1 A ˜ vec ( P k ) and vec 1 M 1 vec ( R k ) , which are required in the PCG method, using the 1, 2, and 3-mode products for tensors and the definition of T. First, from the definitions of A ˜ and X, vec 1 A ˜ vec ( P k ) = vec 1 ( D ¯ X H X D vec ( P k ) ) σ ˜ 2 vec 1 ( X H X vec ( P k ) ) holds. Let D = vec 1 ( diag ( D ) ) , where diag ( D ) returns an m n -dimensional column vector with diagonals of D. Then, D i j k : = λ i ( A ) + λ j ( B ) + λ k ( C ) , where λ i ( A ) , λ j ( B ) , and λ k ( C ) denote the eigenvalues of A , B , and C. Note that ( vec 1 ( D vec ( P k ) ) ) i j k = D i j k ( P k ) i j k since we compute ( D p k ) i = D i i ( p k ) i for i = 1 , 2 , , m n . Using the relation between the Kronecker product and the mode products via vec 1 operator, we compute
vec 1 A ˜ vec ( P ) = D ¯ ( D P k ) × 1 X A H X A + ( D P k ) × 2 X B H X B + ( D P k ) × 3 X C H X C σ ˜ 2 P k × 1 X A H X A + P k × 2 X B H X B + P k × 3 X C H X C ,
where “*” denotes elementwise product.
Next, from the definition of the diagonal matrix M in Equation (5), we easily obtain
M 1 i i = 1 ( D ¯ ) i i ( D ) i i σ ˜ 2 , i = 1 , 2 , , m n .
Here, let M = vec 1 ( diag ( M 1 ) ) . Then it follows that M i j k = 1 / ( D ¯ i j k D i j k σ ˜ 2 ) . vec 1 M 1 vec ( R k ) is computed by
vec 1 M 1 vec ( R k ) = M R k .
As shown in Algorithm 3, the PCG method can be implemented using the preconditioning matrix M and the aforementioned computations, where the linear system A ˜ y ˜ = b ˜ is transformed into A ˜ vec ( Y ˜ ) = vec ( B ˜ ) , where vec ( B ˜ ) : = b ˜ = vec ( Q k × 1 X A H + Q k × 2 X B H + Q k × 3 X C H ) and vec ( Y ˜ ) : = y ˜ . Algorithm 3 requires only small matrices A, B, and C and × m × n tensors X k , R k , P k , and Z k . Therefore the memory requirement is of O ( n 3 ) in the case of n = m = .

4.2. PCG Method by the Schur Decomposition

Firstly, the Schur decomposition of T is T : = Q R Q H , where R and Q are upper triangular and unitary matrices, respectively. Then,
T T T σ ˜ 2 I m n v = q k ( Q R Q H ) H ( Q R Q H ) σ ˜ 2 I m n v = q k R H R σ ˜ 2 I m n ( Q H v ) = Q H q k .
This linear system denotes A ˜ y ˜ = b ˜ , where A ˜ : = R H R σ ˜ 2 I m n , y ˜ : = Q H v , and b ˜ : = Q H q k . The PCG method for A ˜ y ˜ = b ˜ is shown in Algorithm 2. R and Q are obtained from the complex Schur decomposition of A , B , and C as follows: R = I n I m R A + I n R B I + R C I m I and Q = Q C Q B Q A from the definition of T, where A = Q A R A Q A H , B = Q B R B Q B H , and C = Q C R C Q C H by the Schur decomposition of A , B , and C.
Algorithm 3: PCG method over tensor space for the 5-th line of Algorithm 1 [Proposed inner algorithm using the eigendecomposition]
1:  Choose an initial tensor X 0 R × m × n and P 0 = O × m × n , and an initial scalar β 0 = 0 ;
2:  R 0 = ( Q k × 1 X A H + Q k × 2 X B H + Q k × 3 X C H )
    [ D ¯ { ( D X 0 ) × 1 X A H X A + ( D X 0 ) × 2 X B H X B + ( D X 0 ) × 3 X C H X C }
    σ ˜ 2 ( X 0 × 1 X A H X A + X 0 × 2 X B H X B + X 0 × 3 X C H X C ) ] ;
3:  Z 0 = M R 0 ;
4:  for k = 1 , 2 , , until convergence do
5:   P k = Z k 1 + β k 1 P k 1 ;
6:   P ^ k = D ¯ ( D P k ) × 1 X A H X A + ( D P k ) × 2 X B H X B + ( D P k ) × 3 X C H X C
    σ ˜ 2 P k × 1 X A H X A + P k × 2 X B H X B + P k × 3 X C H X C ;
7:   α k = Z k 1 , R k 1 / P k 1 , P ^ k ;
8:   X k = X k 1 + α k P k ;
9:   R k = R k 1 α k P ^ k ;
10:   Z k = M R k 1 ;
11:    β k = Z k , R k / Z k 1 , R k 1 ;
12:  end for
13:  Obtain an approximate solution Y ˜ X k ;
14:  V = Y ˜ × 1 X A + Y ˜ × 2 X B + Y ˜ × 3 X C ;

4.2.1. Preconditioning Matrix

A preconditioning matrix for A ˜ y ˜ = b ˜ satisfies the conditions in Section 4.1.1. Therefore, we propose the preconditioning matrix based on the Schur decomposition
M : = D ¯ R D R σ ˜ 2 I m n ,
where D R is a diagonal matrix with diagonals of R. Since M is also the diagonal matrix, the above second conditions are satisfied. Moreover, if T is symmetric, R is a diagonal matrix, that is, R = D R . Therefore M = A ˜ in the case of the symmetric matrix T. From this, we expect that the preconditioning matrix M is effective if T is not symmetric but almost symmetric.

4.2.2. Efficient Implementation of Algorithm 2 by the Schur Decomposition

We show the computations of vec 1 A ˜ vec ( P k ) and vec 1 M 1 vec ( R k ) for the PCG method over tensor space using the 1, 2, and 3-mode products for tensors and the definition of T. First, from the definitions of A ˜ and R, we have vec 1 A ˜ vec ( P k ) = vec 1 R H R vec ( P k ) σ ˜ 2 vec ( P k ) . Therefore,
vec 1 A ˜ vec ( P k ) = P k × 1 R A H R A + P k × 2 R B H R B + P k × 3 R C H R C σ ˜ 2 P k .
Next, from M = D ¯ R D R σ ˜ 2 I m n , we easily obtain
M 1 i i = 1 ( D ¯ R ) i i ( D R ) i i σ ˜ 2 , i = 1 , 2 , , m n .
Similarly to Section 4.1.2, let D = vec 1 ( diag ( D R ) ) and M = vec 1 ( diag ( M 1 ) ) . Then, we have M i j k = 1 / ( D ¯ i j k D i j k σ ˜ 2 ) , where D i j k = ( R A ) i i + ( R B ) j j + ( R C ) k k . vec 1 ( M 1 vec ( R k ) ) is computed by (7).
As shown in Algorithm 4, the PCG method can be implemented using the preconditioning matrix M and the aforementioned computations, where the linear system A ˜ y ˜ = b ˜ is transformed into A ˜ vec ( Y ˜ ) = vec ( B ˜ ) , where vec ( B ˜ ) : = b ˜ = vec ( Q k × 1 Q A + Q k × 2 Q B + Q k × 3 Q C ) and vec ( Y ˜ ) : = y ˜ . Algorithm 4 just requires small matrices A, B, and C and × m × n tensors X k , R k , P k , and Z k , namely, do not require large matrix T. Therefore the memory requirement is of O ( n 3 ) in the case of n = m = .
Algorithm 4: PCG method over tensor space for the 5-th line of Algorithm 1 [Proposed inner algorithm using the Schur decomposition]
1:  Choose an initial tensor X 0 R × m × n and P 0 = O × m × n , and an initial scalar β 0 = 0 ;
2:  R 0 = ( Q k × 1 Q A + Q k × 2 Q B + Q k × 3 Q C )
    X 0 × 1 R A H R A + X 0 × 2 R B H R B + X 0 × 3 R C H R C σ ˜ 2 X 0 ;
3:  Z 0 = M R 0 ;
4:  for k = 1 , 2 , , until convergence do
5:   P k = Z k 1 + β k 1 P k 1 ;
6:   P ^ k = P k × 1 R A H R A + P k × 2 R B H R B + P k × 3 R C H R C σ ˜ 2 P k ;
7:   α k = Z k 1 , R k 1 / P k 1 , P ^ k ;
8:   X k = X k 1 + α k P k ;
9:   R k = R k 1 α k P ^ k ;
10:   Z k = M R k ;
11:   β k = Z k , R k / Z k 1 , R k 1 ;
12:  end for
13:  Obtain an approximate solution Y ˜ X k ;
14:   V = Y ˜ × 1 Q A + Y ˜ × 2 Q B + Y ˜ × 3 Q C ;

5. Numerical Experiments

This section provides results of numerical experiments using Algorithm 1 with Algorithm 3 and Algorithm 1 with Algorithm 4. There are the two purposes of this experiments: (1) to confirm convergence to the singular value of T nearest to the shift by Algorithm 1, and (2) to confirm the effectiveness of the proposed precondition matrix in Algorithms 3 and 4. For comparison, the results using Algorithms 3 and 4 in the case of M = I are also given as the results by the CG method. All the initial guesses of Algorithms 1, 3, and 4 are tensors with random numbers. The stopping criteria used in Algorithm 1 was β k | | e k T s MAX k | | < 10 8 , where s MAX k is the eigenvector corresponding to the maximum eigenvalue of T ˜ k and e k denotes the k-th canonical basis for k dimensional vector space. Algorithms 3 and 4 were stopped when either the relative residual | | R k | | / | | B ˜ | | < 10 12 or the maximum number of iterations k > 20 , 000 were satisfied.
All computations were carried out using MATLAB R2021a version on a workstation with Xeon processor 3.7 GHz and 128 GB of RAM.
In the following subsection, we show the results computing the 5-th maximum, median, and 5-th minimum singular values σ of the test matrices T. For all the cases, for the first purpose, we set the shift value in Algorithm 1 as
σ ˜ = σ 10 2 ,
where σ ˜ ’s and σ ’s are the perturbed singular values of T and the aforementioned singular values computed by the svd function in MATLAB, respectively.
Test matrices T in Equation (1) are obtained from a seven-point central difference discretization of the PDE (2) in over an ( n + 1 ) × ( n + 1 ) × ( n + 1 ) grid. The test matrices T in Equation (1), whose size is n 3 × n 3 , are generated from
A = B = C , A : = 1 h 2 a M 1 + 1 2 h b M 2 + 1 3 c I n ,
where h = 1 / ( n + 1 ) , M 1 and M 2 are symmetric and skew-symmetric matrices given below.
M 1 = 2 1 1 2 1 1 2 1 1 2 R n × n , M 2 = 0 1 1 0 1 1 0 1 1 0 R n × n .

Numerical Results

In all tables, the number of iterations of the shift-and-invert Lanczos method (“the Lanczos method" hereafter) and the average of the number of iterations of the CG or PCG method based on the eigendecomposition or the Schur decomposition are summarized. “Not converged” denotes Algorithm 3 or 4 did not converge.
We show the first results in the case of almost symmetric matrix with a = c = 1 and b = 0.01 in Equation (9) for the shift (8). From Table 1, Table 2 and Table 3, the numbers of iterations of Lanczos methods using any inner algorithms were almost the same. Focusing on the effectiveness of the proposed preconditioning matrix M, the numbers of iterations of both PCG methods were less than 19 regardless of the size of T. On the other hand, the numbers of iterations of both CG methods linearly increased depending on the size of T. From these facts, the preconditioning matrix M is effective in the case of almost symmetric matrix T. Moreover, the number of iterations of the shift and invert Lanczos method for the median singular value is larger than the number for other singular values since the distance between the maximum and second maximum singular values of ( T T T σ ˜ 2 I m n ) 1 for the median singular value of T is closer than the cases of other singular values.
Here, the running time of Table 1 is summarized in Table 4. All the running time by the PCG method were less than the time by the CG method. Moreover, the running time by the PCG methods of Algorithms 3 and 4 were similar since the computational complexities of these algorithms are similar. Thus, the running time is strongly correlated with the number of iterations of Algorithms 3 and 4.
In addition, convergence histories of n = 15 in Table 1 and Table 2 are shown in Figure 1 and Figure 2. Figure 2 displays that the relative residual norms unsteadily decreased when the number of iterations of the shift and invert Lanczos method is not small.
Next, we show the second results in the case of slightly symmetric matrix with a = c = 1 and b = 0.1 in Equation (9) for the shift (8). From Table 5, both PCG methods did not converge for computing the 5-th maximum singular values of slightly symmetric matrix T. It seems that the linear system for T T T σ ˜ 2 I m n is ill-conditioned since 10 2 in the shift (8) is much less than the 5-th maximum singular values of the matrix. In Appendix A, we show the results using relative shift without the effect of the magnitude of the singular values. Table 6 shows Algorithms 1 and 4, that is, the algorithms based on Schur decomposition, was more robust than Algorithms 1 and 3, that is, the algorithms based on the eigendecomposition. From Table 7, both PCG methods converged regardless of n, and the numbers of iterations of both PCG methods were less than the number of iterations of both CG method. Namely, it seems that the preconditioning matrix M can be effective in the case of a slightly symmetric matrix T when computing the 5-th minimum and median singular values of T.
Finally, we show the third results in the case of marginally symmetric matrix with a = c = 1 and b = 0.2 in Equation (9) for the shift (8). Both PCG methods did not converge for computing the 5-th maximum singular values of T as shown in Table 8, similarly to Table 5. Moreover, computing the median singular values of T sometimes did not converge from Table 9. In Table 10, all methods converged for the 5-th minimum singular value of T. The numbers of iterations by the PCG method with the proposed preconditioning matrix were less than the number of iterations by the CG method in most cases. It seems that the preconditioning matrix M can be effective in the case of the marginally symmetric matrix T when computing the 5-th minimum singular values of T.

6. Conclusions

We considered computing an arbitrary singular value of a tensor sum. The shift-and-invert Lanczos method and the PCG method reconstructed over tensor space. We proposed the preconditioning matrices which are the following two diagonal matrices: (1) whose diagonals of the eigenvalues by the eigendecomposition, and (2) whose diagonals of the upper diagonal matrix by the Schur decomposition. This preconditioning matrix can be effective if the tensor sum is almost symmetric.
From numerical results, we confirmed that the proposed method reduces memory requirements without any conditions. The numbers of iterations of the PCG method by the proposed preconditioning matrix were reduced in most cases of the almost and slightly symmetric matrix. Moreover, the numbers of iterations of the PCG method by the proposed preconditioning matrix were also reduced in certain cases of the marginally symmetric matrix.
For future work, we will consider a robust preconditioning matrix for slightly or marginally symmetric tensor sum, a suitable preconditioning matrix for non-symmetric tensor sum, parallel implementations of the proposed algorithms, and finding real-life applications.

Author Contributions

Conceptualization, A.O. and T.S.; investigation, A.O. and T.S.; software, A.O.; validation, T.S.; writing—original draft, A.O.; writing—review and editing, A.O. and T.S. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by JSPS KAKENHI Grant Numbers: JP21K17754, JP20H00581, JP20K20397, JP17H02829.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
PCGPreconditioned Conjugate Gradient
CGConjugate Gradient
PDEPartial Differential Equation

Appendix A

This appendix gives the numerical results in the case of the 5-th maximum and the median singular values of slightly and marginally symmetric matrices by the relative shift
σ ˜ = σ 10 2 σ ,
where σ ’s are the singular values of T computed by the svd function in MATLAB. The condition of the numerical experiments except for the setting of the shift is the same as the experiments in Section 5.
Firstly, we show the results in the case of slightly symmetric matrix with a = c = 1 and b = 0.1 in Equation (9) for the shift (A1) in Table A1 and Table A2. Computing the 5-th and the median singular values of the slightly symmetric matrix using the shift (A1), the number of iterations of both PCG methods is much less than the number of iterations of both CG methods.
Secondly, Table A3 and Table A4 are the results in the case of marginally symmetric matrix with a = c = 1 and b = 0.2 in Equation (9) for the shift (A1). From Table A3 and Table A4, both PCG methods converged faster than both CG method using the relative shift. Moreover, the PCG method by Algorithm 4 is more robust than the PCG method by Algorithm 3.
Consequently, when we compute the 5-th maximum and the median singular values of the slightly symmetric matrix, the numerical experiments of Section 5 and Appendix A imply that the proposed preconditioning matrix can work in the case of a suitable shift.
Table A1. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (A1).
Table A1. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (A1).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n5541.0517.0535.8515.0
10784.0720.0777.0713.0
154162.0422.04154.0417.0
207223.7723.07197.3712.0
255383.6524.05307.4515.0
306522.7624.06426.5613.0
Table A2. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (A1).
Table A2. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (A1).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n523139.32338.023105.92314.0
10101081.01021.0101074.41022.0
15215201.62121.0213470.42114.0
20177333.1(Not converged.)176242.41716.0
251116,034.11132.01114,360.61114.0
30(Not converged.)1223.0(Not converged.)1217.0
Table A3. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (A1).
Table A3. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (A1).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n5543.0525.0538.4517.0
10790.0729.0779.0713.0
154174.5432.04152.5462.0
207253.0733.07198.1715.0
255403.0535.05319.0519.0
306626.0636.06441.2616.0
Table A4. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (A1).
Table A4. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (A1).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n524138.3(Not converged.)25115.62317.0
10101479.0(Not converged.)101119.21018.0
15214506.62134.0213787.02116.0
20178603.3(Not converged.)176604.51721.0
251118,267.0(Not converged.)1114,991.61130.0
30(Not converged.)1335.0(Not converged.)1331.0

References

  1. Ohashi, A.; Sogabe, T. On computing maximum/minimum singular values of a generalized tensor sum. Electron. Trans. Numer. Anal. 2015, 43, 244–254. [Google Scholar]
  2. Ohashi, A.; Sogabe, T. On computing the minimum singular value of a tensor sum. Special Matrices 2019, 7, 95–106. [Google Scholar] [CrossRef]
  3. Bai, Z.; Demmel, J.; Dongarra, J.; Ruhe, A.; van der Vorst, H.A. Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide; SIAM: Philadelphia, PA, USA, 2000. [Google Scholar]
  4. Beik, F.P.A.; Jbilou, K.; Najafi-Kalyani, M.; Reichel, L. Golub–Kahan bidiagonalization for ill-conditioned tensor equations with applications. Numer. Algorithms 2020, 84, 1535–1563. [Google Scholar] [CrossRef]
  5. Rezaie, M.; Moradzadeh, A.; Kalateh, A.N. Fast 3D inversion of gravity data using solution space priorconditioned lanczos bidiagonalization. J. Appl. Geophys. 2017, 136, 42–50. [Google Scholar] [CrossRef]
  6. Zhong, H.X.; Xu, H. Weighted Golub-Kahan-Lanczos bidiagonalization algorithms. Electron. Trans. Numer. Anal. 2017, 47, 153–178. [Google Scholar] [CrossRef]
  7. Garber, D.; Hazan, E.; Jin, C.; Musco, C.; Netrapalli, P.; Sidford, A. Faster Eigenvector Computation via Shift-and-Invert Preconditioning. In Proceedings of the 33rd International Conference on Machine Learning, New York, NY, USA, 19–24 June 2016; pp. 2626–2634. [Google Scholar]
  8. Huang, W.Q.; Lin, W.W.; Lu, H.H.S.; Yau, S.T. SIRA: Integrated shift—Invert residual Arnoldi method for graph Laplacian matrices from big data. J. Comput. Appl. Math. 2019, 346, 518–531. [Google Scholar] [CrossRef] [Green Version]
  9. Katrutsa, A.; Botchev, M.; Oseledets, I. Practical shift choice in the shift-and-invert Krylov subspace evaluations of the matrix exponential. arXiv 2019, arXiv:1909.13059. [Google Scholar]
  10. Xu, Z. Gradient descent meets shift-and-invert preconditioning for eigenvector computation. Adv. Neural. Inf. Process. Syst. 2018, 31, 2825–2834. [Google Scholar]
  11. Yue, S.F.; Zhang, J.J. An extended shift-invert residual Arnoldi method. Comput. Appl. Math. 2021, 40, 1–15. [Google Scholar] [CrossRef]
  12. Zemaityte, M.; Tisseur, F.; Kannan, R. Filtering Frequencies in a Shift-and-invert Lanczos Algorithm for the Dynamic Analysis of Structures. SIAM J. Sci. Comput. 2019, 41, B601–B624. [Google Scholar] [CrossRef] [Green Version]
  13. Zhong, H.X.; Chen, G.L.; Shen, W.Q. Shift and invert weighted Golub-Kahan-Lanczos bidiagonalization algorithm for linear response eigenproblem. J. Comput. Anal. Appl. 2019, 26, 1169–1178. [Google Scholar]
  14. Van Loan, C.F.; Golub, G.H. Matrix Computations; Johns Hopkins University Press: Baltimore, MD, USA, 1983. [Google Scholar]
  15. Kolda, T.G.; Bader, B.W. Tensor Decompositions and Applications. SIAM. Rev. 2008, 51, 455–500. [Google Scholar] [CrossRef]
Figure 1. Convergence histories with relative residual norm of the Lanczos method for the 5-th max. singular value of the almost symmetric matrix T whose size is n = 15 .
Figure 1. Convergence histories with relative residual norm of the Lanczos method for the 5-th max. singular value of the almost symmetric matrix T whose size is n = 15 .
Axioms 10 00211 g001
Figure 2. Convergence histories with relative residual norm of the Lanczos method for the median singular value of the almost symmetric matrix T whose size is n = 15 .
Figure 2. Convergence histories with relative residual norm of the Lanczos method for the median singular value of the almost symmetric matrix T whose size is n = 15 .
Axioms 10 00211 g002
Table 1. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
Table 1. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n5443.0416.0435.8415.0
10490.0417.0486.5415.0
153134.7317.03128.0317.0
204180.0417.04169.3417.0
253225.7317.03211.0317.0
303273.0317.03252.3317.0
Table 2. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
Table 2. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n515139.31513.015109.11515.0
1071081.0713.07943.7715.0
15415201.63915.0394858.14117.0
20138339.5416.046847.3415.0
25(Not converged.)4817.0(Not converged.)4817.0
30(Not converged.)819.0(Not converged.)815.0
Table 3. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
Table 3. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n53124.7313.0397.9314.6
103748.1313.03654.9314.1
1534872.5314.934531.1316.7
203775.5313.832358.8315.0
2531494.0316.831472.0316.8
3032158.0316.832088.0315.0
Table 4. Running time (seconds) of the Lanczos method using the CG/PCG method in the case of the 5-th max. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
Table 4. Running time (seconds) of the Lanczos method using the CG/PCG method in the case of the 5-th max. singular value of almost symmetric matrix T with a = c = 1 and b = 0.01 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
Lanczos with CGLanczos with PCGLanczos with CGLanczos with PCG
n50.1140.0710.0950.061
100.5780.0950.5660.086
151.4020.1371.4480.192
206.6390.4376.4460.335
2513.6320.55812.6860.432
3035.1211.14533.2030.836
Table 5. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
Table 5. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n5450.0(Not converged.)437.8(Not converged.)
104100.0(Not converged.)487.3(Not converged.)
153152.0(Not converged.)3131.0(Not converged.)
204205.5(Not converged.)4172.0(Not converged.)
253262.0(Not converged.)3230.0(Not converged.)
303306.7(Not converged.)3259.0(Not converged.)
Table 6. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
Table 6. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n513231.8(Not converged.)13193.11373.0
1061582.7655.061272.7629.0
153012,174.6(Not converged.)3111061.93197.0
20413,777.8(Not converged.)48799.3(Not converged.)
25(Not converged.)(Not converged.)(Not converged.)83116.0
30(Not converged.)(Not converged.)(Not converged.)2745.0
Table 7. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
Table 7. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of slightly symmetric matrix T with a = c = 1 and b = 0.1 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n53195.5359.03163.9368.4
103949.8358.03771.2344.0
153616.0363.03597.3394.5
203859.8363.032999.8357.0
2531695.7363.031559.33114.4
3032388.0363.032262.0346.1
Table 8. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
Table 8. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th max. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n5452.8(Not converged.)439.8(Not converged.)
104107.0(Not converged.)493.0(Not converged.)
153162.0(Not converged.)3135.7(Not converged.)
204217.8(Not converged.)4206.3(Not converged.)
253271.3(Not converged.)3248.0(Not converged.)
303334.0(Not converged.)3260.0(Not converged.)
Table 9. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
Table 9. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the median of singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n511481.4(Not converged.)10355.6(Not converged.)
1062123.8(Not converged.)101550.8104262.4
15(Not converged.)(Not converged.)(Not converged.)(Not converged.)
201513,268.7896358.9710019.9108160.0
25(Not converged.)(Not converged.)(Not converged.)(Not converged.)
30(Not converged.)901150.0(Not converged.)2811,764.0
Table 10. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
Table 10. Number of iterations of the Lanczos method and the average of numbers of iterations of the CG/PCG method in the case of the 5-th min. singular value of marginally symmetric matrix T with a = c = 1 and b = 0.2 in (9) for the shift (8).
MethodAlgorithms 1 and 3 (by Eigendecompn.)Algorithms 1 and 4 (by Schur Decompn.)
LanczosCGLanczosPCGLanczosCGLanczosPCG
n55313.95295.85213.951077.8
1031247.83187.031150.633048.8
153650.03201.03606.73177.0
203929.336151.33847.83162.8
2531800.33205.031683.73249.0
3032585.731118.632371.33181.0
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ohashi, A.; Sogabe, T. Numerical Algorithms for Computing an Arbitrary Singular Value of a Tensor Sum. Axioms 2021, 10, 211. https://doi.org/10.3390/axioms10030211

AMA Style

Ohashi A, Sogabe T. Numerical Algorithms for Computing an Arbitrary Singular Value of a Tensor Sum. Axioms. 2021; 10(3):211. https://doi.org/10.3390/axioms10030211

Chicago/Turabian Style

Ohashi, Asuka, and Tomohiro Sogabe. 2021. "Numerical Algorithms for Computing an Arbitrary Singular Value of a Tensor Sum" Axioms 10, no. 3: 211. https://doi.org/10.3390/axioms10030211

APA Style

Ohashi, A., & Sogabe, T. (2021). Numerical Algorithms for Computing an Arbitrary Singular Value of a Tensor Sum. Axioms, 10(3), 211. https://doi.org/10.3390/axioms10030211

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop