The variable projection method is a classical and efficient method for solving separable nonlinear least squares (SNLLS) problems. However, it is hard to handle the constrained SNLLS problems since the explicit form of Jacobian matrix is required in each iteration. In this paper, we propose a secant variable projection (SVP) method, which employs a rank-one update to estimate the Jacobian matrices. The main advantages of our method are efficiency and ease of applicability to constrained SNLLS problems. The local convergence of our SVP method is also analyzed. Finally, some data fitting and image processing problems are solved to compare the performance of our proposed method with the classical variable projection method. Numerical results illustrate the efficiency and stability of our proposed SVP method in solving the SNLLS problems arising from the blind deconvolution problems.

Consider a primal-dual pair of SDP and assume that both primal and dual are either weakly feasible or weakly infeasible. Under such circumstances, the pair might have a nonzero duality gap. Nevertheless, there are arbitrary small perturbations which makes the perturbed primal-dual pair strongly feasible thus zeroing the duality gap. We conduct a perturbation analysis on the pair. We fix two positive definite matrices, \(I_p\) and \(I_d\), say, (typically the identity matrices), and perturb the original primal and dual problem by \(\alpha I_p\) and \(\beta I_d\), respectively to recover interior feasibility, where \(\alpha\) and \(\beta\) are positive numbers. Then we analyze the behavior of the optimal value of the perturbed pair when the perturbation is reduced to zero keeping the proportion of \(\alpha\) and \(\beta\) constant. It will be shown that the optimal value of the perturbed problem converges to a value between the primal and dual optimal values of the original problem. Furthermore, we demonstrate that the infeasible interior-point algorithm for SDP generates a sequence converging to a number between the primal and dual optimal values, even in the presence of a nonzero duality gap.

This is a joint work with Bruno F. Lourenço, Masakazu Muramatsu and Takayuki Okuno.

Tensor decomposition has been widely used for the dimensional reduction and extraction of the meaningful latent features of high dimensional tensor data. In many applications, the underlying data ensemble is nonnegative and consequently the nonnegative tensor decomposition is proposed to achieve additive parts-based representation and to learn more physically interpretable results. As the corresponding tensor optimization problem has computational difficulty due to nonconvex, together with sparse, smooth, graph based Tikhonov regularization, the construction and analysis of the reliable, efficient and robust algorithms are required. Under the framework of block coordinate descent method, we aim to present a new iterative algorithm which is based on the modulus type variable transformation. The theoretical analysis of the proposed method is discussed. Numerical experiments including the synthetic data and image data show the efficiency and superiority of the proposed method comparing with the state-of-the-art methods.

We will study the numerical behaviour of GMRES and block GMRES for coefficient matrices with multiple and clustered eigenvalues and try to analyze it.

In our beating heart simulations based on a fluid-structure coupling finite element model, several factors generate indefinite matrices from Lagrange multiplier approaches dealing with the constraints by the incompressibility of continuums and the contacts between surfaces. In this talk, the strategies in constructing the preconditioner, the iterative method, and the parallelization for our problem are introduced. Furthermore, missing theories for estimating the efficiency of our preconditioner and the convergence of the Krylov subspace iterative method are also discussed.

The preconditioned modified Hermitian/skew-Hermitian splitting (PMHSS) iteration method and the corresponding preconditioning technique can achieve satisfactory results for solving optimal control problems governed by Poisson’s equation. We explore the feasibility of such a method and preconditioner for solving optimization problems constrained by the more complicated Stokes system. Theoretical results demonstrate that the PMHSS iteration method is convergent because the spectral radius of the iterative matrix is less than \(\sqrt{2}/2\). Additionally, the PMHSS preconditioner still clusters eigenvalues on a unitary segment. It guarantees that the convergence of the PMHSS iteration method and preconditioning is independent of not only discretizing mesh size, but also of the Tikhonov regularization parameter. A more effective preconditioner is proposed based on the PMHSS preconditioner. The proposed preconditioner avoids the inner iterations when solving saddle point systems appearing in the generalized residual equations. Furthermore, it is still convergent and maintains its independence of parameter and mesh size.

We consider verifying all eigenvalues in a prescribed interval of a generalized Hermitian eigenproblem. A given eigenproblem can be reduced by a contour integral method into a generalized eigenproblem consisting of complex moments. This kind of method motivated us to develop verification methods that are efficient in parallel. Evaluating all errors in computing the complex moments, we develop verification methods for the eigenproblems. Based on the truncation errors for the trapezoidal rule of numerical quadrature, we derive a numerically computable enclosure of the complex moment. Then, we take into account numerical errors of the quadrature and rigorously enclose the moment matrices. Our error bounds are interval arithmetic-friendly and enables to determine the required number of quadrature points beforehand, according to the degree of precision one needs and the cost one can pay. Each quadrature point gives rise to a linear system, and its structure enables us to develop an efficient technique to verify the approximate solution, which does not need a numerical matrix inverse. Numerical experiments show that the proposed method outperforms a standard method and infer that the proposed method will be potentially efficient in parallel.

This is joint work with Akira Imakura and Akitoshi Takayasu.

In this talk, when \(A\) and \(B\) are \(\{1;1\}\)-quasiseparable matrices, we consider the structured generalized relative eigenvalue condition numbers of the pair \((A, B)\) with respect to relative perturbations of the parameters defining \(A\) and \(B\) in the quasiseparable and the Givens-vector representations of these matrices. A general expression is derived for the condition number of the generalized eigenvalue problems of the pair \((A, B)\), where \(A\) and \(B\) are any differentiable function of a vector of parameters with respect to perturbations of such parameters. Moreover, the explicit expressions of the corresponding structured condition numbers with respect to the quasiseparable and Givens-vector representation via tangents for \(\{1; 1\}\)-quasiseparable matrices are derived. Our proposed condition numbers can be computed efficiently by utilizing the recursive structure of quasiseparable matrices. We investigate relationships between various condition numbers of structured generalized eigenvalue problem when \(A\) and \(B\) are \(\{1;1\}\)-quasiseparable matrices. Numerical results show that there are situations in which the unstructured condition number can be much larger than the structured ones.

In this talk, we study low rank matrix approximation (NLRM) for nonnegative matrices arising from many data mining and pattern recognition applications. Our approach is different from classical nonnegative matrix factorization (NMF) which has been studied for some time. For a given nonnegative matrix, the usual NMF approach is to determine two nonnegative low rank matrices such that the distance between their product and the given nonnegative matrix is as small as possible. However, the proposed NLRM approach is to determine a nonnegative low rank matrix such that the distance between such matrix and the given nonnegative matrix is as small as possible. There are two advantages. (i) The minimized distance can be smaller. (ii) The proposed method can identify important singular basis vectors, while this information may not be obtained in the classical NMF. Numerical results are reported to demonstrate the performance of the proposed method. Several extensions and research works are also presented.

Preconditioners are used to accelerate the convergence of iterative linear solvers. When the linear system \(Ax=b\), changes a little, let's say to \((A+B)x=c\), where \(B\) has low rank, the preconditioner obtained from \(A\) loses effectivity. Instead of computing a new preconditioner from scratch we proposed to update it by bordering the coefficient matrix, and adding new entries to the right-hand side, in such a way that the bordered system is equivalent to the system \((A+B)x=c\). Then the old preconditoner can be reused and updated computing who the added blocks should be incorporated to the preconditioner of the bordered matrix. In addition, the application of the preconditioner should be changed to include the effect of the updating. Moreover, the technique can be used to deal with other situations as rank deficiency of the coefficient matrix, which prevent to obtain a preconditioner using factorization techniques. Also, if the coefficient matrix has a low-rank skew symmetric part, this technique allows to compute a preconditioner for the symmetric part and then include the skew-symmetric in the preconditioning step. In this talk I will review the technique and show some applications to the above mentioned problems.

We apply novel inner-iteration preconditioned Krylov subspace methods to the interior-point (IP) algorithm for linear programming (LP). Inner-iteration preconditioners proposed by Morikuni and Hayami enable us to overcome the severe ill-conditioning of linear equations solved in the final phase of IP iterations. The Krylov subspace methods do not suffer from rank-deficiency and therefore no preprocessing is necessary even if rows of the constraint matrix are not linearly independent. Via these methods, a new IP recurrence is proposed to omit one matrix-vector product per step. Extensive numerical experiments are conducted over diverse instances of 140 LP problems including the Netlib, QAPLIB, Mittelmann and Atomizer Basis Pursuit collections. The largest problem has 434,580 unknowns. It turns out that our implementation is more robust than the standard public domain solvers SeDuMi, SDPT3 and the LSMR iterative solver in PDCO without increasing CPU time. The proposed IP method based on iterative solvers succeeds in solving a fairly large number of LP instances from benchmark libraries under the standard stopping criteria. The work also presents a fairly extensive benchmark test for several renowned solvers including direct and iterative solvers.

We describe computational methods for the solution of \( \ell_p \)-\( \ell_q \) minimization problems of the form \( \min_{x}\{\|Ax-b\|^p_p+\mu\|Lx\|_q^q \} \) with \( 0 \lt q \lt 1 \) and \( 0 \lt p \lt 2 \). Here \( A \) and \( L \) are matrices, \( x \) and \( b \) are vectors, and the parameter \( \mu > 0 \) balances the influence of the two terms. Applications include image restoration and parameter estimation in linear regression. In the latter application, the method is capable of simultaneously producing sparse estimates and dealing with outliers and heavy-tailed error distributions.

Dynamic mode decomposition (DMD) is an efficient method for time series analysis. From a viewpoint of matrix computations, the DMD computes eigenvalues and the corresponding eigenvectors of a given matrix associated with time series data. The basic tools are the singular value decomposition (SVD) and the Rayleigh-Ritz procedure. In this talk, we consider a random matrix due to noisy datasets. In such a case, the DMD computes perturbed eigenpairs including the noise. To reduce the errors, the total least squares method is introduced in the DMD in an implicit form. This algorithm has an excellent performance in terms of the noise cancellation. My recent study focuses on theoretical justification of such kinds of DMD and its generalization. The main result is convergence theory in the probabilistic sense. More specifically, for sample data gradually increasing, it can be proved that the computed eigenpairs contaminated with the random noise converge to the exact ones under a statistical model. In this talk, I will talk about this statistical model and the strong convergence of the computed eigenpairs.

Consider the symmetric singular systems which arise, for instance, in the discretization of partial differential equations for electromagnetic fields using the edge-based finite element method.

When the right-hand side vector is not in the range space of the coefficient matrix, we will compare the convergence of GMRES and Range Restricted GMRES (RRGMRES) methods, which are robust against rounding errors, by numerical experiments and the theory.

We will report the reason that GMRES does not sometimes converge enough unlike the convergence theory about the GMRES method. We will analyze the convergence of the method theoretically, which we have proposed as the solver for the least squares problems when the Hessenberg matrix becomes ill-conditioned, and report the numerical results. Finally, we will show that the range space component of the RRGMRES methods is equivalent to the GMRES methods applied to the nonsingular systems, and will analyze its convergence.

In order to speedup the convergence of GMRES methods for least square problems, different preconditioners including diagonal scaling, robust incomplete factorization, incomplete Givens orthogonalization and Greville's method are studied and analyzed in details, as well as the theoretical properties and implementations. Numerical experiments show they are efficient for the solution of least square problem of suitable size. For large sparse problems, flexible preconditioners updated with the process of the outer iteration are suggested.

For solving large sparse, overdetermined, and inconsistent system of linear equations by iteration methods, by further reconstructing the randomized extended Kaczmarz method proposed by Zouzias and Freris in 2013 (SIAM J. Matrix Anal. Appl. 34(2013), 773–793), we propose a partially randomized extended Kaczmarz method. When the coefficient matrix is assumed to be of full column rank, we prove the convergence and derive an upper bound for the expected convergence rate of the partially randomized extended Kaczmarz method. This bound could be smaller than that of the randomized extended Kaczmarz method under certain conditions. Moreover, with numerical results we show that the partially randomized extended Kaczmarz method can be much more effective than the randomized extended Kaczmarz method.

This talk is concerned with solving some structured multi-linear systems, especially focusing on the equations whose coefficient tensors are M-tensors, or called M-equations for short. We prove that a nonsingular M-equation with a positive right-hand side always has a unique positive solution. Several iterative algorithms are proposed for solving multi-linear nonsingular M-equations, generalizing the classical iterative methods and the Newton method for linear systems. Furthermore, we apply theM-equations to some nonlinear differential equations and the inverse iteration for spectral radii of nonnegative tensors.

Adaptive optics is a commonly used technique to correct the phase distortions caused by the Earth's atmosphere to improve the image quality of the ground-based imaging systems. However, the observed images still suffer from the blur caused by the adaptive optics residual wavefront. We propose a model for reconstructing the residual phase in high-resolution from a sequence of low-resolution deformable mirror data. Our model is based on the turbulence statistics and the Taylor frozen flow hypothesis. A tomography problem for the phase distortions from different altitudes is solved to get a high-quality phase reconstruction. The associated joint optimization problem was solved by an alternating minimization method which results in a high-resolution reconstruction algorithm with adaptive wind velocities. Numerical simulations are carried out to show the effectiveness of our approach.

Joint work with Rihuan Ke (University of Cambridge), Roland Wagner, and Ronny Ramlau (RICAM).

We will reconsider the definition of breakdown of GMRES for singular systems in Brown and Walker (1997). We will also correct the proof of the convergence of GMRES for inconsistent singular systems in Hayami and Sugihara, M. (2011). This is joint work with Dr. Kota Sugihara.

We propose using greedy and randomized Kaczmarz inner-iterations as preconditioners for the right-preconditioned flexible GMRES method to solve consistent linear systems, with a parameter tuning strategy for adjusting the number of inner iterations and the relaxation parameter. We also present theoretical justifications of the right-preconditioned flexible GMRES for solving consistent linear systems. Numerical experiments on overdetermined and underdetermined linear systems show that the proposed method is superior to the GMRES method preconditioned by NE-SOR inner iterations in terms of total CPU time.

Algebraic iterative reconstruction methods — such as ART (Kaczmarz), SART, and SIRT — are popular for CT reconstruction; they adapt to any scanner geometry and handle limited-angle and limited-data problems. The iterative methods utilize projection and backprojection operations in each iteration. We represent the projection by a matrix \(A\) while, in principle, we represent the backprojection by \(B = A^T\) (the transpose of \(A\)).

When the iterative methods are implemented with a focus on computational efficiency, different discretization schemes are used for the forward projection and the backprojection. Hence, there is a mismatch between the backprojection matrix \(B\) and the transposed \(A^T\) of the forward projection matrix. The use of such an unmatched pair has two consequences: the accuracy (compared to when using a matched pair) deteriorates, and the iteration may fail to converge.

In this survey talk, I illustrate these issues with recent theoretical and computational results, and I present a few novel approaches to "fixing" the non-convergence – either by modifying the iterative method or by using the AB- and BA-GMRES methods.

This is joint work with Yiqiu Dong, Tommy Elfving, Ken Hayami, Michiel Hochstenbach, Keiichi Morikuni, and Nicolai Riis.

Randomized low-rank matrix approximation has become an extremely successful approach for efficiently computing a low-rank approximation of matrices. In particular the paper by Halko, Martinsson, and Tropp (SIREV 2011) contains extensive analysis, and has made it a very popular method. The typical complexity for a rank-r approximation of mxn matrices is \(O(mn\log n+(m+n)r^2)\) for dense matrices. In this talk I will first describe a fast and stable algorithm requiring \(O(mn\log n+r^3)\) flops. I also plan to report on recent developments in randomized algorithms: (a subset of) efficient randomized rank estimation, linear solvers, QR factorization, and approximating matrices whose singular values do not decay.