We propose a new strategy for gradient-based unconstrained optimization, involving two parallel sequences of iterates that cooperate to determine their stepsizes via a Twin-Step principle. Rather than minimizing the objective function individually, the algorithm selects steplengths that minimize the Euclidean distance between the two gradient based processes occurring simultaneously at each iteration. The theoretical analysis shows that the convergence of the mutual distance is governed by the angle between the search directions. In particular the effectiveness of the overall process degrades as the directions approach parallelism. To ensure robustness against collinearity, we introduce a hybrid framework, Twin-ABB_min, which switches to the Adaptive Barzilai–Borwein method when the geometric cooperation becomes ineffective. Extensive and very promising numerical results evidence that the Twin phase creates favorable initial conditions for subsequent BB-type iterations.

In the context of low-precision computation for the training of neural networks with the gradient descent method (GD), the occurrence of deterministic rounding errors often leads to stagnation or adversely affects the convergence of the optimizers. We analyze the impact of stochastic floating-point roundoff errors on the convergence behavior of GD with unbiased approximation of gradients. We show that the employment of unbiased stochastic rounding ( SR ) may only yield a similar or slower convergence than exact arithmetic. Additionally, given the assumption of unbiased gradient approximation and SR , the convergence of GD only depends on the variance in evaluating the gradient. When we consider an unbiased gradient approximation, our theoretical claims are also applicable to the stochastic gradient descent method (SGD). We demonstrate, both theoretically and empirically, that the biased stochastic rounding signed-SR_ε can potentially achieve faster convergence when using a smaller stepsize compared to SR with a larger stepsize.

We study Kaczmarz type methods to solve consistent linear matrix equations. We first present a block Kaczmarz (BK) method that employs a deterministic cyclic row selection strategy. Assuming that the associated coefficient matrix has full column or row rank, we derive matrix formulas for a cycle of this BK method. Moreover, we propose a greedy randomized block Kaczmarz (GRBK) method and further extend it to a relaxed variant (RGRBK) and a deterministic counterpart (MWRBK). We establish the convergence properties of the proposed methods. Numerical tests verify the theoretical findings, and we apply the proposed methods to color image restoration problems.

The Kaczmarz method is successfully used for solving discretizations of linear inverse problems, especially in computed tomography where it is known as ART. Practitioners often observe and appreciate its fast convergence in the first few iterations, leading to the same favorable semi-convergence that we observe for simultaneous iterative reconstruction methods. While the latter methods have symmetric and positive definite iteration operators that facilitate their analysis, the operator in Kaczmarz's method is nonsymmetric and it has been an open question so far to understand this fast initial convergence. We perform a spectral analysis of Kaczmarz's method that gives new insight into its (often fast) initial behavior. We also carry out a statistical analysis of how the data noise enters the iteration vectors, which sheds new light on the semi-convergence. Our results are illustrated with several numerical examples.

In Parts I and II of this series of papers, three new methods for the computation of eigenvalues of singular pencils were developed: rank-completing perturbations, rank-projections, and augmentation. It was observed that a straightforward structure-preserving adaption for symmetric pencils was not possible and it was left as an open question how to address this challenge. In this Part III, it is shown how the observed issue can be circumvented by using Hermitian perturbations. This leads to structure-preserving analogues of the three techniques from Parts I and II for Hermitian pencils (including real symmetric pencils) as well as for related structures. It is an important feature of these methods that the sign characteristic of the given pencil is preserved. As an application, it is shown that the resulting methods can be used to solve systems of bivariate polynomials.

When training neural networks with low-precision computation, rounding errors often cause stagnation or are detrimental to the convergence of the optimizers; in this paper we study the influence of rounding errors on the convergence of the gradient descent method for problems satisfying the Polyak-Lojasiewicz inequality. Within this context, we show that, in contrast, biased stochastic rounding errors may be beneficial since choosing a proper rounding strategy eliminates the vanishing gradient problem and forces the rounding bias in a descent direction. Furthermore, we obtain a bound on the convergence rate that is stricter than the one achieved by unbiased stochastic rounding. The theoretical analysis is validated by comparing the performances of various rounding strategies when optimizing several examples using low-precision fixed-point number formats.

A subspace method is introduced to solve large-scale trace ratio problems. This approach is matrix-free, requiring only the action of the two matrices involved in the trace ratio. At each iteration, a smaller trace ratio problem is addressed in the search subspace. Additionally, the algorithm is endowed with a restarting strategy, that ensures the monotonicity of the trace ratio value throughout the iterations. The behavior of the approximate solution is investigated from a theoretical viewpoint, extending existing results on Ritz values and vectors, as the angle between the search subspace and the exact solution approaches zero. Numerical experiments in multigroup classification show that this new subspace method tends to be more efficient than iterative approaches relying on (partial) eigenvalue decompositions at each step.

We present block variants of the discrete empirical interpolation method (DEIM); as a particular application, we will consider a CUR factorization. The block DEIM algorithms are based on the concept of the maximum volume of submatrices and a rank-revealing QR factorization. We also present a version of the block DEIM procedures, which allows for adaptive choice of block size. The results of the experiments indicate that the block DEIM algorithms exhibit comparable accuracy for low-rank matrix approximation compared to the standard DEIM procedure. However, the block DEIM algorithms also demonstrate potential computational advantages, showcasing increased efficiency in terms of computational time.

Recently, three numerical methods for the computation of eigenvalues of singular matrix pencils, based on a rank-completing perturbation, a rank-projection, or an augmentation were developed. We show that all three approaches can be generalized to treat singular polynomial eigenvalue problems. The common denominator of all three approaches is a transformation of a singular into a regular matrix polynomial whose eigenvalues are a disjoint union of the eigenvalues of the singular polynomial, called true eigenvalues, and additional fake eigenvalues. The true eigenvalues can then be separated from the fake eigenvalues using information on the corresponding left and right eigenvectors. We illustrate the approaches on several interesting applications, including bivariate polynomial systems and ZGV points.

The limited memory steepest descent method (LMSD, Fletcher, 2012) for unconstrained optimization problems stores a few past gradients to compute multiple stepsizes at once. We review this method and propose new variants. For strictly convex quadratic objective functions, we study the numerical behavior of different techniques to compute new stepsizes. In particular, we introduce a method to improve the use of harmonic Ritz values. We also show the existence of a secant condition associated with LMSD, where the approximating Hessian is projected onto a low-dimensional space. In the general nonlinear case, we propose two new alternatives to Fletcher’s method: first, the addition of symmetry constraints to the secant condition valid for the quadratic case; second, a perturbation of the last differences between consecutive gradients, to satisfy multiple secant equations simultaneously. We show that Fletcher’s method can also be interpreted from this viewpoint.

We present a new restricted SVD-based CUR (RSVD-CUR) factorization for matrix triplets (A, B, G) that aims to extract meaningful information by providing a low-rank approximation of the three matrices using a subset of their rows and columns. The proposed method utilizes the discrete empirical interpolation method (DEIM) to select the subset of rows and columns from the orthogonal and nonsingular matrices obtained through a restricted singular value decomposition of the matrix triplet. We explore the relationships between a DEIM type RSVD-CUR factorization, a DEIM type CUR factorization, and a DEIM type generalized CUR decomposition, and provide an error analysis that establishes the accuracy of the RSVD-CUR decomposition within a factor of the approximation error of the restricted singular value decomposition of the given matrices. The RSVD-CUR factorization can be used in applications that require approximating one data matrix relative to two other given matrices. We discuss two such applications, namely multi-view dimension reduction and data perturbation problems where a correlated noise matrix is added to the input data matrix. Our numerical experiments demonstrate the advantages of the proposed method over the standard CUR approximation in these scenarios.

Standard multiparameter eigenvalue problems (MEPs) are systems of linear ‐parameter square matrix pencils. Recently, a new form of multiparameter eigenvalue problems has emerged: a rectangular MEP (RMEP) with only one multivariate rectangular matrix pencil, where we are looking for combinations of the parameters for which the rank of the pencil is not full. Applications include finding the optimal least squares autoregressive moving average (ARMA) model and the optimal least squares realization of autonomous linear time‐invariant (LTI) dynamical system. For linear and polynomial RMEPs, we give the number of solutions and show how these problems can be solved numerically by a transformation into a standard MEP. For the transformation we provide new linearizations for quadratic multivariate matrix polynomials with a specific structure of monomials and consider mixed systems of rectangular and square multivariate matrix polynomials. This numerical approach seems computationally considerably more attractive than the block Macaulay method, the only other currently available numerical method for polynomial RMEPs.

A CUR factorization is often utilized as a substitute for the singular value decomposition (SVD), especially when a concrete interpretation of the singular vectors is challenging. Moreover, if the original data matrix possesses properties like nonnegativity and sparsity, a CUR decomposition can better preserve them compared to the SVD. An essential aspect of this approach is the methodology used for selecting a subset of columns and rows from the original matrix. This study investigates the effectiveness of one-round sampling and iterative subselection techniques and introduces new iterative subselection strategies based on iterative SVDs. One provably appropriate technique for index selection in constructing a CUR factorization is the discrete empirical interpolation method (DEIM). Our contribution aims to improve the approximation quality of the DEIM scheme by iteratively invoking it in several rounds, in the sense that we select subsequent columns and rows based on the previously selected ones. Thus, we modify A after each iteration by removing the information that has been captured by the previously selected columns and rows. We also discuss how iterative procedures for computing a few singular vectors of large data matrices can be integrated with the new iterative subselection strategies. We present the results of numerical experiments, providing a comparison of one-round sampling and iterative subselection techniques, and demonstrating the improved approximation quality associated with using the latter.

Given an approximate eigenvector, its (standard) Rayleigh quotient and harmonic Rayleigh quotient are two well-known approximations of the corresponding eigenvalue. We propose a new type of Rayleigh quotient, the homogeneous Rayleigh quotient, and analyze its sensitivity with respect to perturbations in the eigenvector. Furthermore, we study the inverse of this homogeneous Rayleigh quotient as stepsize for the gradient method for unconstrained optimization. The notion and basic properties are also extended to the generalized eigenvalue problem.

We study the use of inverse harmonic Rayleigh quotients with target for the stepsize selection in gradient methods for nonlinear unconstrained optimization problems. This not only provides an elegant and flexible framework to parametrize and reinterpret existing stepsize schemes, but it also gives inspiration for new flexible and tunable families of steplengths. In particular, we analyze and extend the adaptive Barzilai–Borwein method to a new family of stepsizes. While this family exploits negative values for the target, we also consider positive targets. We present a convergence analysis for quadratic problems extending results by Dai and Liao (IMA J Numer Anal 22(1):1–10, 2002), and carry out experiments outlining the potential of the approaches.

Generalized eigenvalue problems involving a singular pencil may be very challenging to solve, both with respect to accuracy and efficiency. While Part I presented a rank-completing addition to a singular pencil, we now develop two alternative methods. The first technique is based on a projection onto subspaces with dimension equal to the normal rank of the pencil while the second approach exploits an augmented matrix pencil. The projection approach seems to be the most attractive version for generic singular pencils because of its efficiency, while the augmented pencil approach may be suitable for applications where a linear system with the augmented pencil can be solved efficiently.

When implementing the gradient descent method in low precision, the employment of stochastic rounding schemes helps to prevent stagnation of convergence caused by the vanishing gradient effect. Unbiased stochastic rounding yields zero bias by preserving small updates with probabilities proportional to their relative magnitudes. This study provides a theoretical explanation for the stagnation of the gradient descent method in low-precision computation. Additionally, we propose two new stochastic rounding schemes that trade the zero bias property with a larger probability to preserve small gradients. Our methods yield a constant rounding bias that, on average, lies in a descent direction. For convex problems, we prove that the proposed rounding methods typically have a beneficial effect on the convergence rate of gradient descent. We validate our theoretical analysis by comparing the performances of various rounding schemes when optimizing a multinomial logistic regression model and when training a simple neural network with an 8-bit floating-point format.

We compare two different linear dimensionality reduction strategies for the multigroup classification problem: the trace ratio method and Fisher's discriminant analysis. Recently, trace ratio optimization has gained in popularity due to its computational efficiency, as well as the occasionally better classification results. However, a statistical understanding is still incomplete. We study and compare the properties of the two methods. Then, we propose a robust version of the trace ratio method, to handle the presence of outliers in the data. We reinterpret an asymptotic perturbation bound for the solution to the trace ratio, in a contamination setting. Finally, we compare the performance of the trace ratio method and Fisher's discriminant analysis on both synthetic and real datasets, using classical and robust estimators.

The discrete empirical interpolation method (DEIM) may be used as an index selection strategy for formulating a CUR factorization. A notable drawback of the original DEIM algorithm is that the number of column or row indices that can be selected is limited to the number of input singular vectors. We propose a new variant of DEIM, which we call L-DEIM, a combination of the strength of deterministic leverage scores and DEIM. This method allows for the selection of a number of indices greater than the number of input singular vectors. Since DEIM requires singular vectors as input matrices, L-DEIM is particularly attractive for example in big data problems when computing a rank-k SVD approximation is expensive even for moderately small k since it uses a lower-rank SVD approximation instead of the full rank-k SVD. We empirically demonstrate the performance of L-DEIM, which despite its efficiency, may achieve comparable results to the original DEIM and even better approximations than some state-of-the-art methods.