Extended Krylov subspace methods for solving Sylvester and Stein tensor equations

Abdeslem Hafid Bentbib; Smahane El-Halouy; El Mostafa Sadek

doi:10.3934/dcdss.2021026

Article Contents

2022, Volume 15, Issue 1: 41-56. Doi: 10.3934/dcdss.2021026

This issue Previous Article A new coupled complex boundary method (CCBM) for an inverse obstacle problem Next Article An efficient D-N alternating algorithm for solving an inverse problem for Helmholtz equation

Extended Krylov subspace methods for solving Sylvester and Stein tensor equations

1.
Laboratory LAMAI, University Cadi Ayyad, Marrakesh, Morocco
2.
Laboratory LabSIPE, ENSA d'El Jadida, University Chouaib Doukkali, El Jadida, Morocco

^* Corresponding author: Smahane El-Halouy
^* Corresponding author: Smahane El-Halouy

Received: July 31, 2020

Revised: December 31, 2020

Early access: March 2021

Published: January 2022

Abstract Full Text(HTML) Figure(3) / Table(8) Related Papers Cited by

Abstract

This paper deals with Sylvester and Stein tensor equations with low rank right hand sides. It proposes extended Krylov-like methods for solving Sylvester and Stein tensor equations. The expressions of residual norms are presented. To show the performance of the proposed approaches, some numerical experiments are given.

Keywords:

Mathematics Subject Classification: 2020 MSC.

Citation:

Full Text(HTML)

Figure 1. Residual norm of the global and the extended global Arnoldi algorithms for example $ 3 $ with the first set of matrices and $ n = 225 $

Download: Full-size image PowerPoint slide

Figure 2. Residual norm of the global and the extended global Arnoldi algorithms for example $ 3 $ with the second set of matrices and $ n = 225 $

Download: Full-size image PowerPoint slide

Figure 3. Residual norm of the extended block and global Arnoldi algorithms for example $ 4 $ with $ n = 300 $

Download: Full-size image PowerPoint slide

Algorithm 1 Extended Block Arnoldi (EBA)

1: Input: $ n\times n $ matrix $ A $, $ n\times r $ matrix $ B $ and an integer $ m $.
2: Compute the QR decomposition of $ \left[ B,A^{-1}B\right]=U_1\Lambda, \mathbb{U}_0=\left[ \;\right] $ for $ j=1,\ldots ,m $
3: Set $ U_j^{(1)} $: first $ r $ columns of $ U_j $, $ U_j^{(2)} $: second $ r $ columns of $ U_j $
4: $ \mathbb{U}_j=\left[ \mathbb{U}_{j-1},U_j\right] $; $ W=\left[ AU_j^{(1)},A^{-1}U_j^{(2)}\right] $
5: Orthogonalize $ W $ w.r. to $ \mathbb{U}_j $ to get $ U_{j+1} $ i.e.,
6: for $ i=1,\ldots ,j $ $ H_{i,j}=U_i^TW $; $ W=W-U_iH_{i,j} $
7: Compute the QR decomposition of $ W=U_{j+1}H_{j+1,j} $

| Show Table

DownLoad: CSV

Algorithm 2 Extended Global Arnoldi (EGA)

1: Input: $ n\times n $ matrix $ A $, $ n\times r $ matrix $ B $ and an integer $ m $.
2: Compute the global QR decomposition of $ \left[ B,A^{-1}B\right]=V_1(\Omega\otimes I_r) $; $ \mathbb{V}_0=\left[ \;\right] $ for $ j=1,\ldots ,m $
3: Set $ V_j^{(1)} $: first $ r $ columns of $ V_j $, $ V_j^{(2)} $: second $ r $ columns of $ V_j $
4: $ \mathbb{V}_j=\left[ \mathbb{V}_{j-1},V_j\right] $; $ U=\left[ AV_j^{(1)},A^{-1}V_j^{(2)}\right] $
5: Orthogonalize $ U $ w.r. to $ \mathbb{V}_j $ to get $ V_{j+1} $ i.e.,
6: for $ i=1,\ldots ,j $ $ H_{i,j}=V_i^T\lozenge U $; $ U=U-V_iH_{i,j} $
7: Compute the QR decomposition of $ U=V_{j+1}(H_{j+1,j}\otimes I_r) $

| Show Table

DownLoad: CSV

Algorithm 3 Sylvester Stein EBA Process

1: Coefficient matrices $ A^{(i)} $, $ i=1,\ldots ,N $, and the right hand side in low rank representation $ \mathcal{B}=\left[B^{(1)},\ldots ,B^{(N)} \right] $.
2: Output: Approximate solutions, $ \mathcal{X}_m $, to the equations (1) / (2).
3: Choose a tolerance $ \epsilon>0 $, a maximum number of iterations $ itermax_i $, a step-size parameter $ k $, set for $ i=1,\ldots ,N $, $ m_i=k $.
4: For $ i=1,\ldots ,N $ do
For $ m_i=k $, construct the orthonormal basis $ \mathbb{V}_{m_i} $ and the restriction matrices $ \mathbb{T}_{m_i} $ by EBA algorithm (1).
5: Compute $ \mathcal{Y}_m $ / $ \mathcal{Z}_m $ the solution of the low dimensional equation (12) / (13).
6: Compute the residual norms $ r_m=\|\mathcal{R}_m\| $ as in theorem (4.1)/ $ q_m=\|\mathcal{Q}_m\| $ as in theorem (4.2).
7: If $ r_m\geq\epsilon $ / $ q_m\geq\epsilon $, set for $ i=1,\ldots ,N $, $ m_i=mi+k $ and go to step 2.
8: The approximate solutions are given by $ \mathcal{X}_m=\mathcal{Y}_m\times_1\mathbb{V}_{m_1}\ldots\times_N\mathbb{V}_{m_N} $ / $ \mathcal{X}_m=\mathcal{Z}_m\times_1\mathbb{V}_{m_1}\ldots\times_N\mathbb{V}_{m_N} $.

| Show Table

DownLoad: CSV

Algorithm 4 Sylvester / Stein EGA Process

1: Input: Coefficient matrices $ A^{(i)} $, $ i=1,\ldots ,N $, and the right hand side in low rank representation $ \mathcal{B}=\left[B^{(1)},\ldots ,B^{(N)} \right] $.
2: Output: Approximate solutions, $ \mathcal{X}_m $, to the equations (1) / (2).
3: Choose a tolerance $ \epsilon>0 $, a maximum number of iterations $ itermax_i $, a step-size parameter $ k $, set for $ i=1,\ldots ,N $, $ m_i=k $.
4: For $ i=1,\ldots ,N $ do
For $ m_i=k $, construct the orthonormal basis $ \mathbb{W}_{m_i} $ and the restriction matrices $ T_{m_i} $ by EGA algorithm (2).
5: Compute $ \mathcal{Y}_m $ / $ \mathcal{Z}_m $ the solutions of the low dimensional equations (16) / (19).
6: Compute the upper bound for the residual norms $ r_m=\left( {\sum_{i=1}^N\|\mathcal{Y}_m\times_iT_{m_i+1,m_i}E_{m_i}^T\|^2}\right)^{1/2} $ as in theorem (4.4)/ $ q_m=\left( \beta_1^2+\sum_{i=1}^N{\beta_2^{(i)}}^2+\sum_{i=1}^N{\beta_3^{(i)}}^2\right)^{1/2} $ as in theorem (4.5).
7: If $ r_m\geq\epsilon $ / $ q_m\geq\epsilon $, set for $ i=1,\ldots ,N $, $ m_i=mi+k $ and go to step 2.
8: The approximate solutions are given by $ \mathcal{X}_m=(\mathcal{Y}_m\otimes\mathcal{I}_R)\times_1\mathbb{W}_{m_1}\ldots\times_N\mathbb{W}_{m_N} $ / $ \mathcal{X}_m=(\mathcal{Z}_m\otimes\mathcal{I}_R)\times_1\mathbb{W}_{m_1}\ldots\times_N\mathbb{W}_{m_N} $.

| Show Table

DownLoad: CSV

Table 1. Example $ 1 $

Algorithm	$ R $	$ n $	Cycles	$ \mathcal{R}_m $	CPU Times (Second)
Algorithm 3	$ 3 $	$ 500 $	$ 5 $	$ 9.34\cdot 10^{-7} $	$ 16.1 $
	$ 5 $	$ 300 $	$ 4 $	$ 4.03\cdot 10^{-7} $	$ 31.38 $
Algorithm 4	$ 3 $	$ 500 $	$ 9 $	$ 6.69\cdot 10^{-7} $	$ 4.57 $
	$ 5 $	$ 300 $	$ 8 $	$ 4.34\cdot 10^{-7} $	$ 3.76 $

| Show Table

DownLoad: CSV

Table 2. Example 2

Algorithm		$ R $	$ n $	Cycles	$ \\|\mathcal{R}_m\\| $	$ \\|\mathcal{X}^{\ast}-\mathcal{X}_m\\| $	CPU Times (Second)
Algorithm 3	(case $ 1 $)	$ 3 $	$ 300 $	$ 2 $	$ 4.08\cdot 10^{-9} $	$ 3.85\cdot 10^{-6} $	$ 1.58 $
	(case $ 2 $)	$ 3 $	$ 300 $	$ 2 $	$ 1.26\cdot 10^{-8} $	$ 2.82\cdot 10^{-7} $	$ 1.53 $
Algorithm 4	(case $ 1 $)	$ 3 $	$ 300 $	$ 6 $	$ 5.82\cdot 10^{-8} $	$ 4.21\cdot 10^{-6} $	$ 1.46 $
	(case $ 2 $)	$ 3 $	$ 300 $	$ 7 $	$ 5.39\cdot 10^{-8} $	$ 4.17\cdot 10^{-8} $	$ 2.22 $
Algorithm 4		$ 5 $	$ 500 $	$ 10 $	$ 4.54\cdot 10^{-7} $	$ - $	$ 26.01 $

| Show Table

DownLoad: CSV

Table 3. Example 3, $ N = 3 $ and $ R = 5 $

Matrices	Algorithm	$ \\|\mathcal{R}_m\\| $	CPU Times (Second)
fdm	Algorithm $ 2 $ [4]	$ 1.45\cdot 10^{-6} $	$ 56.0 $
	Algorithm 4	$ 3.6\cdot 10^{-9} $	$ 0.74 $
poisson	Algorithm $ 2 $ [4]	$ 2.82\cdot 10^{-6} $	$ 14.5 $
	Algorithm 4	$ 2.24\cdot 10^{-9} $	$ 0.67 $

| Show Table

DownLoad: CSV

Table 4. Example $ 4 $

Algorithm	$ R $	$ n $	Cycles	$ \mathcal{R}_m $	CPU Times (Second)
Algorithm 3	$ 3 $	$ 300 $	$ 6 $	$ 1.58\cdot 10^{-6} $	$ 27.49 $
Algorithm 4	$ 3 $	$ 300 $	$ 4 $	$ 1.56\cdot 10^{-7} $	$ 0.83 $
	$ 5 $	$ 300 $	$ 4 $	$ 8.77\cdot 10^{-7} $	$ 0.85 $

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	J. Ballani and L. Grasedyck, A projection method to solve linear systems in tensor format, Numerical Linear Algebra with Applications, 20 (2013), 27–43. doi: 10.1002/nla.1818.
[2]	F. P. A. Beik, F. S. Movahed and S. Ahmadi-Asl, On the Krylov subspace methods based on tensor format for positive definite Sylvester tensor equations, Numerical Linear Algebra with Applications, 23 (2016), 444–466. doi: 10.1002/nla.2033.
[3]	A. H. Bentbib, K. Jbilou and El M. Sadek, On some extended block Krylov based methods for large scale nonsymmetric Stein matrix equations, Mathematics, 5 (2017), 21. doi: 10.3390/math5020021.
[4]	A. H. Bentbib, S. El-Halouy and El M. Sadek, Krylov subspace projection method for Sylvester tensor equation with low rank right-hand side, Numerical Algorithms, (2020), 1–20. doi: 10.1007/s11075-020-00874-0.
[5]	M. Chen and D. Kressner, Recursive blocked algorithms for linear systems with Kronecker product structure, Numerical Algorithms, (2019), 1–18. doi: 10.1007/s11075-019-00797-5.
[6]	D. Calvetti and L. Reichel, Application of ADI iterative methods to the restoration of noisy images, SIAM Journal on Matrix Analysis and Applications, 17 (1996), 165–186. doi: 10.1137/S0895479894273687.
[7]	Z. Chen and L. Lu, A projection method and Kronecker product pre-conditioner for solving Sylvester tensor equations, Science China Mathematics, 55 (2012), 1281-1292. doi: 10.1007/s11425-012-4363-5.
[8]	Z. Chen and L. Lu, A gradient based iterative solutions for Sylvester tensor equations, Mathematical Problems in Engineering, (2013). doi: 10.1155/2013/819479.
[9]	F. Ding and T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Transactions on Automatic Control, 50 (2005), 1216–1221. doi: 10.1109/TAC.2005.852558.
[10]	M. Heyouni and K. Jbilou, An extended block Arnoldi algorithm for large-scale solutions of the continuous-time algebraic Riccati equation, Electronic Transactions on Numerical Analysis, 33 (2009), 53–62.
[11]	M. Heyouni, Extended Arnoldi methods for large low-rank Sylvester matrix equations, Applied Numerical Mathematics, 60 (2010), 1171–1182. doi: 10.1016/j.apnum.2010.07.005.
[12]	T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review, 51 (2009), 455–500. doi: 10.1137/07070111X.
[13]	D. Kressner and C. Tobler, Krylov subspace methods for linear systems with tensor product structure, SIAM Journal on Matrix Analysis and Applications, 31 (2010), 1688–1714. doi: 10.1137/090756843.
[14]	A. N. Langville and W. J. Stewart, A Kronecker product approximate preconditioner for SANs, Numerical Linear Algebra with Applications, 11 (2004), 723–752. doi: 10.1002/nla.344.
[15]	N. Lee and A. Cichocki, Fundamental tensor operations for large-scale data analysis using tensor network formats, Multidimensional Systems and Signal Processing, 29 (2018), 921–960. doi: 10.1007/s11045-017-0481-0.
[16]	T. Li, Q. W. Wang and X. F. Duan, Numerical algorithms for solving discrete Lyapunov tensor equation, Journal of Computational and Applied Mathematics, 370 (2020), 112676. doi: 10.1016/j.cam.2019.112676.
[17]	T. Penzl, et al., A Matlab toolbox for large Lyapunov and Riccati equations, model reduction problems, and linear quadratic optimal control problems, https://www.tu-chemnitz.de/sfb393/lyapack/, (2000).
[18]	T. Stykel and V. Simoncini, Krylov subspace methods for projected Lyapunov equations, Applied Numerical Mathematics, 62 (2012), 35-50. doi: 10.1016/j.apnum.2011.09.007.
[19]	X. Xu and W. Q. Wang, Extending BiCG and BiCR methods to solve the Stein tensor equation, Computers & Mathematics with Applications, 77 (2019), 3117–3127. doi: 10.1016/j.camwa.2019.01.024.