Fractional Tikhonov regularization for linear discrete ill-posed problems

Tikhonov regularization is one of the most popular methods for solving linear systems of equations or linear least-squares problems with a severely ill-conditioned matrix A. This method replaces the given problem by a penalized least-squares problem. The present paper discusses measuring the residual error (discrepancy) in Tikhonov regularization with a seminorm that uses a fractional power of the Moore-Penrose pseudoinverse of AA ^T as weighting matrix. Properties of this regularization method are discussed. Numerical examples illustrate that the proposed scheme for a suitable fractional power may give approximate solutions of higher quality than standard Tikhonov regularization.     

A pivoting strategy for symmetric tridiagonal matrices

The LBL^T factorization of Bunch for solving linear systems involving a symmetric indefinite tridiagonal matrix T is a stable, efficient method. It computes a unit lower triangular matrix L and a block 1 × 1 and 2 × 2 matrix B such that T=LBL^T. Choosing the pivot size requires knowing a priori the largest element σ of T in magnitude. In some applications, it is required to factor T as it is formed without necessarily knowing σ. In this paper, we present a modification of the Bunch algorithm that can satisfy this requirement. We demonstrate that this modification exhibits the same bound on the growth factor as the Bunch algorithm and is likewise normwise backward stable. Copyright © 2005 John Wiley & Sons, Ltd.     

Sign consistent linear programming problems

This article studies linear programming problems in which all minors of maximal order of the coefficient matrix have the same sign. We analyse the relationship between a special structure of the non-degenerate dual feasible bases of a linear programming problem and the structure of its associated matrix. In the particular case in which the matrix has all minors of each order k with the same strict sign ? _k , we provide a dual simplex revised method with good stability properties. In particular, this method can be applied to the totally positive linear programming problems, of great interest in many applications.     

Numerical solution to low rank perturbed Lyapunov equations by the sign function method

We investigate the numerical solution to a low rank perturbed Lyapunov equation A^TX + XA = W via the sign function method (SFM). The sign function method has been proposed to solve Lyapunov equations, see e.g. [1], but here we focus on a framework where the matrix A has a special structure, i. e. A = B + UCV^T , where B is a blockdiagonal matrix and UCV^T is a low rank perturbation. We show that this structure can be kept throughout the sign function iteration but the rank of the perturbation doubles per iteration. Therefore, we apply a low rank approximation to the perturbation in order to keep its numerical rank small. We compare the standard SFM with its structure preserving variant presented in this paper by means of numerical examples from viscously damped mechanical systems. (© 2016 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Chapter 4:algebraic sets,polynomials, and other functions

Let T be a linear operator on the vector space V ofn×n matrices over a field F. We discuss two types of problems in this chapter. First, what can we say about T if we assume that T maps a given algebraic set such as the special linear group into itself? Second, let p(x) be a polynomial function (such as det) on V into F. What can we say about T if Tpreserves p(x), i.e. p(T(X)) = p(X) for all X in V?     

Hyper-Sparsity in the Revised Simplex Method and How to Exploit it

The revised simplex method is often the method of choice when solving large scale sparse linear programming problems, particularly when a family of closely-related problems is to be solved. Each iteration of the revised simplex method requires the solution of two linear systems and a matrix vector product. For a significant number of practical problems the result of one or more of these operations is usually sparse, a property we call hyper-sparsity. Analysis of the commonly-used techniques for implementing each step of the revised simplex method shows them to be inefficient when hyper-sparsity is present. Techniques to exploit hyper-sparsity are developed and their performance is compared with the standard techniques. For the subset of our test problems that exhibits hyper-sparsity, the average speedup in solution time is 5.2 when these techniques are used. For this problem set our implementation of the revised simplex method which exploits hyper-sparsity is shown to be competitive with the leading commercial solver and significantly faster than the leading public-domain solver.     

Using successive approximations for improving the convergence of GMRES method

In this paper, our attention is concentrated on the GMRES method for the solution of the system (I–T)x=b of linear algebraic equations with a nonsymmetric matrix. We perform m pre-iterations y _l+1 =T _yl +b before starting GMRES and put y _m for the initial approximation in GMRES. We derive an upper estimate for the norm of the error vector in dependence on the mth powers of eigenvalues of the matrix T Further we study under what eigenvalues lay-out this upper estimate is the best one. The estimate shows and numerical experiments verify that it is advisable to perform pre-iterations before starting GMRES as they require fewer arithmetic operations than GMRES. Towards the end of the paper we present a numerical experiment for a system obtained by the finite difference approximation of convection-diffusion equations.     

On Copositive Programming and Standard Quadratic Optimization Problems

A standard quadratic problem consists of finding global maximizers of a quadratic form over the standard simplex. In this paper, the usual semidefinite programming relaxation is strengthened by replacing the cone of positive semidefinite matrices by the cone of completely positive matrices (the positive semidefinite matrices which allow a factorization FF ^T where F is some non-negative matrix). The dual of this cone is the cone of copositive matrices (i.e., those matrices which yield a non-negative quadratic form on the positive orthant). This conic formulation allows us to employ primal-dual affine-scaling directions. Furthermore, these approaches are combined with an evolutionary dynamics algorithm which generates primal-feasible paths along which the objective is monotonically improved until a local solution is reached. In particular, the primal-dual affine scaling directions are used to escape from local maxima encountered during the evolutionary dynamics phase.     

Worst case behavior of the steepest edge simplex method

At each nondegenerate iteration of the steepest-edge simplex method one moves from a vertex of the polyhedron, P, of feasible points to an adjacent vertex along an edge that is steepest with respect to the linear objective function ψ. In this paper we show how to construct a sequence of linear programs (P_n,ψ_n) in n variables for which the number of iterations required by the steepest edge simplex method is 2ⁿ−1.     

Comparison of subdominant eigenvalues of some linear search schemes

The subdominant eigenvalue of the transition probability matrix of a Markov chain is a determining factor in the speed of transition of the chain to a stationary state. However, these eigenvalues can be difficult to estimate in a theoretical sense. In this paper we revisit the problem of dynamically organizing a linear list. Items in the list are selected with certain unknown probabilities and then returned to the list according to one of two schemes: the move-to-front scheme or the transposition scheme. The eigenvalues of the transition probability matrix Q of the former scheme are well-known but those of the latter T are not. Nevertheless the transposition scheme gives rise to a reversible Markov chain. This enables us to employ a generalized Rayleigh-Ritz theorem to show that the subdominant eigenvalue of T is at least as large as the subdominant eigenvalue of Q.     

Solution of sparse rectangular systems using LSQR and CRAIG

We examine two iterative methods for solving rectangular systems of linear equations: LSQR for over-determined systemsAx b, and Craig's method for under-determined systemsAx = b. By including regularization, we extend Craig's method to incompatible systems, and observe that it solves the same damped least-squares problems as LSQR. The methods may therefore be compared on rectangular systems of arbitrary shape.Various methods for symmetric and unsymmetric systems are reviewed to illustrate the parallels. We see that the extension of Craig's method closes a gap in existing theory. However, LSQR is more economical on regularized problems and appears to be more reliable if the residual is not small.In passing, we analyze a scaled augmented system associated with regularized problems. A bound on the condition number suggests a promising direct method for sparse equations and least-squares problems, based on indefiniteLDL ^T factors of the augmented matrix.Dedicated to Professor Åke Björck in honor of his 60th birthdayPresented at the 12th Householder Symposium on Numerical Algebra, Lake Arrowhead, California, June 1993.Partially supported by Department of Energy grant DE-FG03-92ER25117, National Science Foundation grant DMI-9204208, and Office of Naval Research grant N00014-90-J-1242.     

An LDLT factorization based ADI algorithm for solving large-scale differential matrix equations

Large-scale differential matrix equations appear in many applications like optimal control of partial differential equations, balanced truncation model order reduction of linear time varying systems etc. Here, we will focus on matrix Riccati differential equations (RDE). Solving such matrix valued ordinary differential equations (ODE) is a highly storage and time consuming process. Therefore, it is necessary to develop efficient solution strategies minimizing both. We present an LDL^T factorization based ADI method for solving algebraic Lyapunov equations (ALE) arising in the innermost iteration during the application of Rosenbrock ODE solvers to RDEs. We show that the LDL^T-type decomposition avoids complex arithmetic, as well as cancellation effects arising from indefinite right hand sides of the ALEs appearing in the classic ZZ^T based approach. Additionally, a certain number of linear system solves can be saved within the ADI algorithm by reducing the number of column blocks in the right hand sides while the full accuracy of the standard low-rank ADI is preserved. (© 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Least-squares solutions of multi-valued linear operator equations in Hilbert spaces

Let M be a linear manifold in H₁ H₂, where H₁, and H₂ are Hilbert spaces. Two notions of least-squares solutions for the multi-valued linear operator equation (inclusion) y ε M(x) are introduced and investigated. The main results include (i) equivalent conditions for least-squares solvability, (ii) properties of a least-squares solution, (iii) characterizations of the set of all least-squares solutions in terms of algebraic operator parts and generalized inverses of linear manifolds, and (iv) best approximation properties of generalized inverses and operator parts of multi-valued linear operators. The principal tools in this investigation are an abstract adjoint theory, orthogonal operator parts, and orthogonal generalized inverses of linear manifolds in Hilbert spaces.     

Preconditioned HSS method for large multilevel block Toeplitz linear systems via the notion of matrix‐valued symbol

We perform a spectral analysis of the preconditioned Hermitian/skew‐Hermitian splitting (PHSS) method applied to multilevel block Toeplitz linear systems in which the coefficient matrix T_n(f) is associated with a Lebesgue integrable matrix‐valued function f. When the preconditioner is chosen as a Hermitian positive definite multilevel block Toeplitz matrix T_n(g), the resulting sequence of PHSS iteration matrices M_n belongs to the generalized locally Toeplitz class. In this case, we are able to compute the symbol ?(f,g) describing the asymptotic eigenvalue distribution of M_nwhen n→∞ and the matrix size diverges. By minimizing the infinity norm of the spectral radius of the symbol ?(f,g), we are also able to identify effective PHSS preconditioners T_n(g) for the matrix T_n(f). A number of numerical experiments are presented and commented, showing that the theoretical results are confirmed and that the spectral analysis leads to efficient PHSS methods. Copyright © 2015 John Wiley & Sons, Ltd.     

SQP algorithms for solving Toeplitz matrix approximation problem

Given an n × n matrix F, we find the nearest symmetric positive semi‐definite Toeplitz matrix T to F. The problem is formulated as a non‐linear minimization problem with positive semi‐definite Toeplitz matrix as constraints. Then a computational framework is given. An algorithm with rapid convergence is obtained by l₁ Sequential Quadratic Programming (SQP) method. Copyright © 2002 John Wiley & Sons, Ltd.     

Solution of sparse quasi-square rectangular systems by Gaussian elimination

We present a general method for the linear least-squares solutionof overdetermined and underdetermined systems. The method isparticularly efficient when the coefficient matrix is quasi-square,that is when the number of rows and number of columns is almostthe same. The numerical methods for linear least-squares problemsand minimum-norm solutions do not generally take account ofthis special characteristic. The proposed method is based onLU factorization of the original quasi-square matrix A, assumingthat A has full rank. In the overdetermined case, the LU factorsare used to compute a basis for the null space of A^T. The right-handside vector b is then projected onto this subspace and the least-squaressolution is obtained from the solution of this reduced problem.In the case of underdetermined systems, the desired solutionis again obtained through the solution of a reduced system.The use of this method may lead to important savings in computationaltime for both dense and sparse matrices. It is also shown inthe paper that, even in cases where the matrices are quite small,sparse solvers perform better than dense solvers. Some practicalexamples that illustrate the use of the method are included.     

Continuous-time linear programming problems revisited: A perturbation approach

The class of continuous-time linear programming problems under the assumption that the constraints are satisfied almost everywhere in the time interval [0,?T]?is taken into account in this article. Under this assumption, its corresponding discretized problems cannot be formulated by equally dividing the time interval [0,?T]?as subintervals of [0,?T]?. In this article, we also introduce the perturbed continuous-time linear programming problems to prove the strong duality theorem when the constraints are assumed to be satisfied a.e. in [0,?T]?.     

On constrained generalized inverses of matrices and their properties

A matrix G is called a generalized inverse (g-invserse) of matrix A if AGA = A and is denoted by G = A ⁻. Constrained g-inverses of A are defined through some matrix expressions like E(AE)⁻, (FA)⁻ F and E(FAE)⁻ F. In this paper, we derive a variety of properties of these constrained g-inverses by making use of the matrix rank method. As applications, we give some results on g-inverses of block matrices, and weighted least-squares estimators for the general linear model.     

Tensor-on-Tensor Regression

I propose a framework for the linear prediction of a multiway array (i.e., a tensor) from another multiway array of arbitrary dimension, using the contracted tensor product. This framework generalizes several existing approaches, including methods to predict a scalar outcome from a tensor, a matrix from a matrix, or a tensor from a scalar. I describe an approach that exploits the multiway structure of both the predictors and the outcomes by restricting the coefficients to have reduced PARAFAC/CANDECOMP rank. I propose a general and efficient algorithm for penalized least-squares estimation, which allows for a ridge (L₂) penalty on the coefficients. The objective is shown to give the mode of a Bayesian posterior, which motivates a Gibbs sampling algorithm for inference. I illustrate the approach with an application to facial image data. An R package is available at https://github.com/lockEF/MultiwayRegression.     

Optimal Hierarchies for Quadrilateral Surfaces

Multiresolution representation of quadrilateral surface approximation (MRQSA) is a useful representation for progressive graphics transmission in networks. Based on two requirements: (1) minimum mean square error and (2) fixed reduction ratio between levels, this paper first transforms the MRQSA problem into the problem of solving a sequence of near-Toeplitz tridiagonal linear systems. Employing the matrix perturbation technique, the MRQSA problem can be solved using about 24mn floating-point operations, i.e. linear time, if we are given a polygonal surface with (2m–1)×(2n–1) points. A numerical stability analysis is also given. To the best of our knowledge, this is the first time that such a linear algebra approach has been used for solving the MRQSA problem. Some experimental results are carried out to demonstrate the applicability of the proposed method.