Search

The present paper is concerned with developing tensor iterative Krylov subspace methods to solve large multi-linear tensor equations. We use the T-product for two tensors to define tensor tubal global Arnoldi and tensor tubal global Golub-Kahan bidiagonalization algorithms. Furthermore, we illustrate how tensor-based global approaches can be exploited to solve ill-posed problems arising from recovering blurry multichannel (color) images and videos, using the so-called Tikhonov regularization technique, to provide computable approximate regularized solutions. We also review a generalized cross-validation and discrepancy principle type of criterion for the selection of the regularization parameter in the Tikhonov regularization. Applications to image sequence processing are given to demonstrate the efficiency of the algorithms.

Patch-based low-rank approximation (PLRA) via truncated singular value decomposition is a powerful and effective tool for recovering the underlying low-rank structure in images. Generally, it first performs an approximate nearest neighbors (ANN) search algorithm to group similar patches into a collection of matrices with reshaping them as vectors. The inherent correlation among similar patches makes these matrices have a low-rank structure. Then the singular value decomposition (SVD) is used to derive a low-rank approximation of each matrix by truncating small singular values. However, the conventional implementation of patch-based low-rank image restoration suffers from high computational cost of the ANN search and full SVD. To address this limitation, we propose a fast approximation method that accelerates the computation of PLRA using multiple kd-trees and Lanczos approximation. The basic idea of this method is to exploit an index kd-tree built from patch samples of the observed image and several small kd-trees built from overlapping regions of the image to accelerate the search for similar patches, and apply the Lanczos bidiagonalization procedure to obtain a fast low-rank approximation of patch matrix without computing the full SVD. Experimental results on image denoising and inpainting tasks demonstrate the efficiency and accuracy of our method.

The magnetotelluric (MT) technique is an electromagnetic geophysical method, which is widely used as a complementary to seismic surveys for exploration of hydrocarbon reservoirs. In the inversion process, the method of matrix inverse calculation has a considerable effect on the speed of the inversion and the quality of obtained models. Lanczos Bidiagonalization (LB) method has been reported to be a fast and efficient approach for solving the inversion problems. In this study, we employ LB method for inverting large-scale 2D MT data. In LB algorithm, the full set of equations is replaced by a dimensionally reduced system of equations. As a result, the speed of the solution procedure is increased, while the original problem is solved with a high accuracy. In addition, we employ active constraint balancing approach for determining the optimum regularization parameter. The advantage of the method is that for highly resolvable parameters, a small value of the Lagrangian multiplier is assigned, and vice versa. The results of the synthetic data inversion show that both methods require equal computer memory but LB method is faster and more reliable than conjugate gradient method. The proposed approach is also applied to inverse real MT data collected from the Kashan area. The Kashan area is the most interesting area for oil and gas exploration of the Central Iran Basin. The inversion results obtained by LB are in a good agreement with the geological structure of the study area and the drilling data.

Randomized methods can be competitive for the solution of problems with a large matrix of low rank. They also have been applied successfully to the solution of large-scale linear discrete ill-posed problems by Tikhonov regularization (Xiang and Zou in Inverse Probl 29:085008, 2013). This entails the computation of an approximation of a partial singular value decomposition of a large matrixAthat is of numerical low rank. The present paper compares a randomized method to a Krylov subspace method based on Golub–Kahan bidiagonalization with respect to accuracy and computing time and discusses characteristics of linear discrete ill-posed problems that make them well suited for solution by a randomized method.

The joint bidiagonalization(JBD) process is a useful algorithm for the computation of the generalized singular value decomposition(GSVD) of a matrix pair. However, it always suffers from rounding errors, which causes the Lanczos vectors to loss their mutual orthogonality. In order to maintain some level of orthongonality, we present a semiorthogonalization strategy. Our rounding error analysis shows that the JBD process with the semiorthogonalization strategy can ensure that the convergence of the computed quantities is not affected by rounding errors and the final accuracy is high enough. Based on the semiorthogonalization strategy, we develop the joint bidiagonalization process with partial reorthogonalization(JBDPRO). In the JBDPRO algorithm, reorthogonalizations occur only when necessary, which saves a big amount of reorthogonalization work compared with the full reorthogonalization strategy. Numerical experiments illustrate our theory and algorithm., 26 pages, 7 figures. arXiv admin note: substantial text overlap with arXiv:1912.08505

Structural model updating (SMU) by sensitivity analysis is a reliable and successful approach to adjusting finite element models of real structures by modal data. However, some challenging issues, including the incompleteness of measured modal parameters, the presence of noise, and the ill-posedness of the inverse problem of SMU may cause inaccurate and unreliable updating results. By deriving a hybrid sensitivity function based on a combination of modal strain and kinetic energies, this article proposes a sensitivity-based SMU method for simultaneously updating the elemental mass and stiffness matrices. The key idea of the proposed sensitivity function originates from the fundamental concepts of the equality-constrained optimization problem and Lagrange multiplier method. By combining the Lanczos bidiagonalization algorithm and Tikhonov regularization technique, a novel hybrid regularization method is presented to solve the ill-posed inverse problem of SMU in a robust manner. This scheme exploits new generalized cross-validation functions directly derived from the outputs of the Lanczos bidiagonalization algorithm to automatically specify the number of iterations and determine an optimal regularization parameter. The accuracy and effectiveness of the presented approaches are verified numerically and experimentally by a steel truss and the full-scale I-40 Bridge. All results show that the proposed methods are effective tools for SMU under incomplete noisy modal data.

Based on the joint bidiagonalization (JBD) process of the matrix pair { A , L } , an iterative regularization algorithm, called JBDQR, is proposed and developed for large scale linear discrete ill-posed problems in general-form Tikhonov regularization. It is proved that the JBDQR iterates take the form of attractive filtered generalized singular value decomposition (GSVD) expansions, where the filters are given explicitly and insightful. This result and a detailed analysis on it show that JBDQR must have the desired semi-convergence property, where the iteration number k plays the role of the regularization parameter. Embedded with the L-curve criterion or the discrepancy principle that is used to estimate the optimal k ⁎ at which the semi-convergence occurs, JBDQR can compute a satisfying good regularized solution. JBDQR is theoretically solid and effective, and it is simple to implement. Numerical experiments confirm our results and the robustness of JBDQR.

This paper is concerned with the solution of severely ill-conditioned linear tensor equations. These kinds of equations may arise when discretizing partial differential equations in many space-dimensions by finite difference or spectral methods. The deblurring of color images is another application. We describe the tensor Golub–Kahan bidiagonalization (GKB) algorithm and apply it in conjunction with Tikhonov regularization. The conditioning of the Stein tensor equation is examined. These results suggest how the tensor GKB process can be used to solve general linear tensor equations. Computed examples illustrate the feasibility of the proposed algorithm.

For large-scale problems, how to establish an algorithm with high accuracy and stability is particularly important. In this paper, the Householder bidiagonalization total least squares (HBITLS) algorithm and nonlinear iterative partial least squares for total least squares (NIPALS-TLS) algorithm were established, by which the same approximate TLS solutions was obtained. In addition, the propagation of the roundoff error for the process of the HBITLS algorithm was analyzed, and the mixed forward-backward stability of these two algorithms was proved. Furthermore, an upper bound of roundoff error was derived, which presents a more detailed and clearer approximation of the computed solution.

We consider the iterative method LSQR for solving $\min_x |Ax-b|_2$. LSQR is based on the Golub--Kahan bidiagonalization process and at every step produces an iterate that minimizes the norm of the...

In this paper we analyze the behavior of the LSQR algorithm for the solution of compact operator equations in Hilbert spaces. We present results concerning existence of Krylov solutions and the rate of convergence in terms of an l p sequence where p depends on the summability of the singular values of the operator. Under stronger regularity requirements we also consider the decay of the error. Finally we study the approximation of the dominant singular values of the operator attainable with the bidiagonal matrices generated by the Lanczos bidiagonalization and the arising low rank approximations. Some numerical experiments on classical test problems are presented.

The periodicity of an object is one of its most important visual characteristics. Recently, several low-rank/sparse matrix decomposition techniques have indicated that a relationship exists between the frequency components of the motion matrix and its decomposition components. This relationship was mostly identified based on empirical evidence without proper analysis, which led to an unclear understanding and poor utilization. This paper attempts to establish the relationship between the periodic components in the motion matrix and its singular value decomposition (SVD) components. The transformation of the periodic components of the motion matrix through QR factorization and Golub–Kahan bidiagonalization, which are the two essential steps of SVD, is thoroughly discussed and analyzed. Furthermore, the two cases of fully and partially periodic motion matrices are considered where the relationships between their frequency components and left singular vectors are established. The performance of the proposed SVD-based scheme is validated using real videos with ground truth. The results show that the proposed scheme correctly estimates the singular vectors that contain the frequency information.

A singular-value-decomposition (SVD) processor having adjustable precision for $8 \times 8$ multiple-input multiple-output precoding systems supporting up to 256-quadrature amplitude modulation (QAM) is designed and implemented. The memory-based architecture that consists of eight processing elements, each having two coordinate rotation digital computer modules, is employed. Golub–Reinsch SVD (GR-SVD) algorithm with Rayleigh quotient shift is used. Thus, two-phase operations are performed including bidiagonalization and implicit shifted QR for diagonalization. The split, deflation, and shift techniques of GR-SVD are helpful to accelerate the diagonalization and reduce the computation complexity. To cover the precision requirements for compact 256-QAM constellation and spatially correlated channels, hybrid datapath representations are used. The threshold for split and deflation can be adjusted and thus the precision of the SVD processor is variable according to the requirements. Although the high precision results in a large gate count, this SVD processor in 40-nm CMOS technology can complete the decomposition of an $8 \times 8$ matrix in 313 clock cycles averagely and is able to provide a throughput rate of 591K matrix/s with good power efficiency.

The quaternion signal model can naturally represent vector-field data in a physically consistent approach, plausibly preserving and exploiting its vector relationships. The Singular Spectrum Analysis (SSA) filter can be extended to the case of multicomponent seismic data via the non-commutative quaternion algebra. Adoption of quaternions requires a shift from real or complex signal processing to quaternion-based signal processing. In this paper, the quaternion-defined seismic data is transformed to the frequency–space ( f − x ) domain via the quaternion Fourier transform (QFT), and rank-reduction can be performed via the quaternion SVD (QSVD) following the SSA methodology. Efficient rank-reduction alternatives, such as Lanczos bidiagonalization via convolutions, are used to allow significant gains in processing time. Simultaneous denoising and reconstruction of synthetic multicomponent seismic records are used to illustrate and compare the quaternionic versions of the SSA against its scalar version.

The exponential growth of big data places a great burden on current computing environment. However, there exists a vast gap in the approaches that can securely and efficiently process the large scale heterogeneous data. This paper, on the basis of transparent computing paradigm, presents a unified approach that coordinates the transparent servers and transparent clients to decompose tensor, a mathematical model widely used in data intensive applications, to a core tensor multiplied with a number of truncated orthogonal bases. The structured, semi-structured as well as structured data are transformed to low-order sub-tensors, which are then encrypted using the Paillier homomorphic encryption scheme on the transparent clients. The cipher sub-tensors are transported to the transparent servers for carrying out the integration and decomposition operations. Three secure decomposition algorithms, namely secure bidiagonalization algorithm, secure singular value decomposition algorithm, and secure mode product algorithm, are presented to generate the bidiagonal matrices, truncated orthogonal bases, and core tensor respectively. The homomorphic operations of the three algorithms are carried out on the transparent servers, while the non-homomorphic operations, namely division and square root, are performed on the transparent clients. Experimental results indicate that the proposed method is promising for secure tensor decomposition for big data.

In semiconductor manufacturing process, lithography simulation modeling is known as an ill-posed problem. A normal solution of the problem is generally insignificant due to measurement constraints. In order to alleviate the difficulty, we introduced a regularization method using a preconditioning technique which consists of scaling and uniformization based on prior information. By regularizing the solution to prior knowledge, an accurate model can be achieved because the solution using truncated singular value decomposition from a few larger singular values becomes a reasonable solution based on the physically appropriate prior knowledge. The augmented implicitly restarted Lanczos bidiagonalization (AIRLB) algorithm is suitable for the purpose of the truncated singular value decomposition from a few larger singular values. Thus, the AIRLB algorithm is important for obtaining the solution in lithography simulation modeling. In this paper, we propose techniques for improving the AIRLB algorithm for the truncated singular value decomposition of large matrices. Specifically, we implement the improvement of the AIRLB algorithm by Ishida et al. Furthermore, instead of using the QR algorithm, we use the orthogonal-qd-with-shift algorithm for the singular value decomposition of the inner small matrix. Several numerical experiments demonstrate that, compared with AIRLB using the original QR algorithm, the proposed improvements provide highly accurate truncated singular value decomposition. For precise discussion, both large-scale sparse matrices and large-scale dense matrices are included in the experiments.

Kinematic relationships between degrees of freedom, also named multi-point constraints, are frequently used in structural mechanics. In this paper, the Craig variant of the Golub-Kahan bidiagonalization algorithm is used as an iterative method to solve the arising linear system with a saddle point structure. The condition number of the preconditioned operator is shown to be close to unity and independent of the mesh size. This property is proved theoretically and illustrated on a sequence of test problems of increasing complexity, including concrete structures enforced with pretension cables and the coupled finite element model of a reactor containment building. The Golub-Kahan algorithm converges in only a small number of steps for all considered test problems and discretization sizes. Furthermore, it is robust in practical cases that are otherwise considered to be difficult for iterative solvers.

The efficacy of the tomographic process depends upon the image reconstruction. Utmost mathematical problems encounter in tomography are systems of large linear equations. Krylov solvers for linear systems have sophisticated and straightforward formulae for the residual norm. Two Krylov solvers CGLS and LSQR are the variations of steep descent. The steep descent is one of the fundamental iterative technique used exclusively for the solution of large sparse square matrices. However, CGLS and LSQR the variations of steep descent also solve least square problems. This work involves the comparison of CGLS and LSQR. CGLS and LSQR are mathematically equivalent, but LSQR is robust and difficult to apply. Large sparse linear least square problem solved by LSQR, that is Krylov space solver in-fact based on Lanczos's bidiagonalization. This work applies said variations (CGLS and LSQR) of the steep descent on a tomographic test problem and compare the two algorithms on the basis of accuracy using MATLAB.

We describe LNLQ for solving the least-norm problem min $\|x\|$ subject to $Ax=b$, using the Golub--Kahan bidiagonalization of $[b\ A]$. Craig's method is known to be equivalent to applying the con...

When solving ill-posed inverse problems, a good choice of the prior is critical for the computation of a reasonable solution. A common approach is to include a Gaussian prior, which is defined by a mean vector and a symmetric and positive definite covariance matrix, and to use iterative projection methods to solve the corresponding regularized problem. However, a main challenge for many of these iterative methods is that the prior covariance matrix must be known and fixed (up to a constant) before starting the solution process. In this paper, we develop hybrid projection methods for inverse problems with mixed Gaussian priors where the prior covariance matrix is a convex combination of matrices and the mixing parameter and the regularization parameter do not need to be known in advance. Such scenarios may arise when data is used to generate a sample prior covariance matrix (e.g., in data assimilation) or when different priors are needed to capture different qualities of the solution. The proposed hybrid methods are based on a mixed Golub–Kahan process, which is an extension of the generalized Golub–Kahan bidiagonalization, and a distinctive feature of the proposed approach is that both the regularization parameter and the weighting parameter for the covariance matrix can be estimated automatically during the iterative process. Furthermore, for problems where training data are available, various data-driven covariance matrices (including those based on learned covariance kernels) can be easily incorporated. Numerical examples from tomographic reconstruction demonstrate the potential for these methods.

In many color image processing and recognition applications, one of the most important targets is to compute the optimal low-rank approximations to color images, which can be reconstructed with a small number of dominant singular value decomposition (SVD) triplets of quaternion matrices. All existing methods are designed to compute all SVD triplets of quaternion matrices at first and then to select the necessary dominant ones for reconstruction. This way costs quite a lot of operational flops and CPU times to compute many superfluous SVD triplets. In this paper, we propose a Lanczos-based method of computing partial (several dominant) SVD triplets of the large-scale quaternion matrices. The partial bidiagonalization of large-scale quaternion matrices is derived by using the Lanczos iteration, and the reorthogonalization and thick-restart techniques are also utilized in the implementation. An algorithm is presented to compute the partial quaternion singular value decomposition. Numerical examples, including principal component analysis, color face recognition, video compression and color image completion, illustrate that the performance of the developed Lanczos-based method for low-rank quaternion approximation is better than that of the state-of-the-art methods.

This work discusses four algorithms for the solution of linear discrete ill-posed problems with several right-hand side vectors. These algorithms can be applied, for instance, to multi-channel image restoration when the image degradation model is described by a linear system of equations with multiple right-hand sides that are contaminated by errors. Two of the algorithms are block generalizations of the standard Golub–Kahan bidiagonalization method with the block size equal to the number of channels. One algorithm uses standard Golub–Kahan bidiagonalization without restarts for all right-hand sides. These schemes are compared to standard Golub–Kahan bidiagonalization applied to each right-hand side independently. Tikhonov regularization is used to avoid severe error propagation. Numerical examples illustrate the performance of these algorithms. Applications include the restoration of color images.

We introduce a new family of saddle-point minimum residual methods for iteratively solving saddle-point systems using a minimum or quasi-minimum residual approach. No symmetry assumptions are made. The basic mechanism underlying the method is a novel simultaneous bidiagonalization procedure that yields a simplified saddle-point matrix on a projected Krylov-like subspace and allows for a monotonic short-recurrence iterative scheme. We develop a few variants, demonstrate the advantages of our approach, derive optimality conditions, and discuss connections to existing methods. Numerical experiments illustrate the merits of this new family of methods.

An iterative method, LSMB, is given for solving $\min_x \|Ax-b\|_2$. LSMB is based on the Golub--Kahan bidiagonalization process and is constructed so that an objective function closely related to ...

We consider the solution of large linear systems of equations that arise from the discretization of ill-posed problems. The matrix has a Kronecker product structure and the right-hand side is contaminated by measurement error. Problems of this kind arise, for instance, from the discretization of Fredholm integral equations of the first kind in two space-dimensions with a separable kernel and in image restoration problems. Regularization methods, such as Tikhonov regularization, have to be employed to reduce the propagation of the error in the right-hand side into the computed solution. We investigate the use of the global GolubKahan bidiagonalization method to reduce the given large problem to a small one. The small problem is solved by employing Tikhonov regularization. A regularization parameter determines the amount of regularization. The connection between global GolubKahan bidiagonalization and Gauss-type quadrature rules is exploited to inexpensively compute bounds that are useful for determining the regularization parameter by the discrepancy principle.

We present a scalability study of Golub-Kahan bidiagonalization for the parallel iterative solution of symmetric indefinite linear systems with a 2x2 block structure. The algorithms have been implemented within the parallel numerical library PETSc. Since a nested inner-outer iteration strategy may be necessary, we investigate different choices for the inner solvers, including parallel sparse direct and multigrid accelerated iterative methods. We show the strong and weak scalability of the Golub-Kahan bidiagonalization based iterative method when applied to a two-dimensional Poiseuille flow and to two- and three-dimensional Stokes test problems.

Summary When inverting potential-field geophysical data an appropriate discretization of the model is required to accurately construct complicated geometries of the causative sources. Rectangular prisms (structured meshes) have limitations to recover and preserve the edges of realistic geological sources. We use an isoparametric finite-element methodology to design an unstructured mesh for use in inverse modeling of gravity data. The calculation of the sensitivity kernel of the forward operator uses Gauss-Legendre quadrature rather than the analytic formulation. For the sake of instability of inversion operator in gravity data and to solve the Tikhonov norms due to improvements in the regularization principle associated with the stabilizing term, the Lanczos bidiagonalization technique is used. This method is accompanied with weighted generalized cross-validation (WGCV) for selecting the optimum amount of the regularization parameter. The depth weighting function is also incorporated in the formulation of the objective function to suppress the impact of shallow features and recover sources at an appropriate depth. Our algorithm is applied to real case study where a gravity survey is used for iron exploration in Yazd province, central Iran. The real example recovers geologically reasonable complex structures.

Focusing inversion of potential field data for the recovery of sparse subsurface structures from surface measurement data on a uniform grid is discussed. For the uniform grid the model sensitivity matrices exhibit block Toeplitz Toeplitz block structure, by blocks for each depth layer of the subsurface. Then, through embedding in circulant matrices, all forward operations with the sensitivity matrix, or its transpose, are realized using the fast two dimensional Fourier transform. Simulations demonstrate that this fast inversion algorithm can be implemented on standard desktop computers with sufficient memory for storage of volumes up to size $n \approx 1M$. The linear systems of equations arising in the focusing inversion algorithm are solved using either Golub Kahan bidiagonalization or randomized singular value decomposition algorithms in which all matrix operations with the sensitivity matrix are implemented using the fast Fourier transform. These two algorithms are contrasted for efficiency for large-scale problems with respect to the sizes of the projected subspaces adopted for the solutions of the linear systems. The presented results confirm earlier studies that the randomized algorithms are to be preferred for the inversion of gravity data, and that it is sufficient to use projected spaces of size approximately $m/8$, for data sets of size $m$. In contrast, the Golub Kahan bidiagonalization leads to more efficient implementations for the inversion of magnetic data sets, and it is again sufficient to use projected spaces of size approximately $m/8$. Moreover, it is sufficient to use projected spaces of size $m/20$ when $m$ is large, $m \approx 50000$, to reconstruct volumes with $n \approx 1M$. Simulations support the presented conclusions and are verified on the inversion of a practical magnetic data set that is obtained over the Wuskwatim Lake region in Manitoba, Canada., Comment: 39 pages. 16 figures

We present an iterative method based on a generalization of the Golub-Kahan bidiagonalization for solving indefinite matrices with a 2 \(\times \) 2 block structure. We focus in particular on our recent implementation of the algorithm using the parallel numerical library PETSc. Since the algorithm is a nested solver, we investigate different choices for parallel inner solvers and show its strong scalability for two Stokes test problems. The algorithm is found to be scalable for large sparse problems.