For a datastream, the change over a short interval is often of low rank. For high throughput information arranged in matrix format, recomputing an optimal SVD approximation after each step is typically prohibitive. Instead, incremental and truncated updating strategies are used, which may not scale for large truncation ranks. Therefore, we propose a set of efficient new algorithms that update a bidiagonal factorization, and which are similarly accurate as the SVD methods. In particular, we develop a compact Householder-type algorithm that decouples a sparse part from a low-rank update and has about half the memory requirements of standard bidiagonalization methods. A second algorithm based on Givens rotations has only about 10 flops per rotation and scales quadratically with the problem size, compared to a typical cubic scaling. The algorithm is therefore effective for processing high-throughput updates, as we demonstrate in tracking large subspaces of recommendation systems and networks, and when compared to well known software such as LAPACK or the incremental SVD.

Linear solvers are major computational bottlenecks in a wide range of decision support and optimization computations. The challenges become even more pronounced on heterogeneous hardware, where traditional sparse numerical linear algebra methods are often inefficient. For example, methods for solving ill-conditioned linear systems have relied on conditional branching, which degrades performance on hardware accelerators such as graphical processing units (GPUs). To improve the efficiency of solving ill-conditioned systems, our computational strategy separates computations that are efficient on GPUs from those that need to run on traditional central processing units (CPUs). Our strategy maximizes the reuse of expensive CPU computations. Iterative methods, which thus far have not been broadly used for ill-conditioned linear systems, play an important role in our approach. In particular, we extend ideas from [1] to implement iterative refinement using inexact LU factors and flexible generalized minimal residual (FGMRES), with the aim of efficient performance on GPUs. We focus on solutions that are effective within broader application contexts, and discuss how early performance tests could be improved to be more predictive of the performance in a realistic environment

We present a GPU implementation of Algorithm NCL, an augmented Lagrangian method for solving large-scale and degenerate nonlinear programs. Although interior-point methods and sequential quadratic programming are widely used for solving nonlinear programs, the augmented Lagrangian method is known to offer superior robustness against constraint degeneracies and can rapidly detect infeasibility. We introduce several enhancements to Algorithm NCL, including fusion of the inner and outer loops and use of extrapolation steps, which improve both efficiency and convergence stability. Further, NCL has the key advantage of being well-suited for GPU architectures because of the regularity of the KKT systems provided by quadratic penalty terms. In particular, the NCL subproblem formulation allows the KKT systems to be naturally expressed as either stabilized or condensed KKT systems, whereas the interior-point approach requires aggressive reformulations or relaxations to make it suitable for GPUs. Both systems can be efficiently solved on GPUs using sparse \ldlt factorization with static pivoting, as implemented in NVIDIA cuDSS. Building on these advantages, we examine the KKT systems arising from NCL subproblems. We present an optimized GPU implementation of Algorithm NCL by leveraging MadNLP as an interior-point subproblem solver and utilizing the stabilized and condensed formulations of the KKT systems for computing Newton steps. Numerical experiments on various large-scale and degenerate NLPs, including optimal power flow, COPS benchmarks, and security-constrained optimal power flow, demonstrate that MadNCL operates efficiently on GPUs while effectively managing problem degeneracy, including MPCC constraints.

Algorithm NCL was devised to solve a class of large nonlinearly constrained optimization problems whose constraints do not satisfy LICQ at a solution. It is mathematically equivalent to the augmented Lagrangian algorithm LANCELOT, which solves a short sequence of bound-constrained subproblems BC_k and has no LICQ difficulties. NCL’s equivalent subproblems NC_k are much bigger and must be solved by a nonlinear interior method (needing first and second derivatives). We study the KKT-type systems arising within nonlinear interior methods when they are applied to the NC_k subproblems. We find that the KKT systems can sometimes be reduced to smaller SQD systems (symmetric quasi-definite) and sometimes to even smaller SPD systems (symmetric positive definite). The smaller systems have proved suitable for GPU implementation within the interior solver MadNLP when it is used by MadNCL to implement Algorithm NCL.

We introduce an iterative solver named MINARES for symmetric linear systems $Ax \approx b$, where $A$ is possibly singular. MINARES is based on the symmetric Lanczos process, like MINRES and MINRES-QLP, but it minimizes $\|Ar_k\|$ in each Krylov subspace rather than $\|r_k\|$, where $r_k$ is the current residual vector. When $A$ is symmetric, MINARES minimizes the same quantity $\|Ar_k\|$ as LSMR, but in more relevant Krylov subspaces, and it requires only one matrix-vector product $Av$ per iteration, whereas LSMR would need two. Our numerical experiments with MINRES-QLP and LSMR show that MINARES is a pertinent alternative on consistent symmetric systems and the most suitable Krylov method for inconsistent symmetric systems. We derive properties of MINARES from an equivalent solver named CAR that is to MINARES as CR is to MINRES, is not based on the Lanczos process, and minimizes $\|Ar_k\|$ in the same Krylov subspace as MINARES. We establish that MINARES and CAR generate monotonic $\|x_k - x_{\star}\|$, $\|x_k - x_{\star}\|_A$ and $\|r_k\|$ when $A$ is positive definite.

For linear systems $Ax=b$ we develop iterative algorithms based on a sketch-and-project approach. By using judicious choices for the sketch, such as the history of residuals, we develop weighting strategies that enable short recursive formulas. The proposed algorithms have a low memory footprint and iteration complexity compared to regular sketch-and-project methods. In a set of numerical experiments the new methods compare well to GMRES, SYMMLQ and state-of-the-art randomized solvers.

Augmented Lagrangian (AL) methods are a well known class of algorithms for solving constrained optimization problems. They have been extended to the solution of saddle-point systems of linear equations. We study an AL (SPAL) algorithm for unsymmetric saddle-point systems and derive convergence and semi-convergence properties, even when the system is singular. At each step, our SPAL requires the exact solution of a linear system of the same size but with an SPD (2,2) block. To improve efficiency, we introduce an inexact SPAL algorithm. We establish its convergence properties under reasonable assumptions. Specifically, we use a gradient method, known as the Barzilai-Borwein (BB) method, to solve the linear system at each iteration. We call the result the augmented Lagrangian BB (SPALBB) algorithm and study its convergence. Numerical experiments on test problems from Navier-Stokes equations and coupled Stokes-Darcy flow show that SPALBB is more robust and efficient than BICGSTAB and GMRES. SPALBB often requires the least CPU time, especially on large systems.

The conjugate gradient (CG) method is a classic Krylov subspace method for solving symmetric positive definite linear systems. We introduce an analogous semi-conjugate gradient (SCG) method for unsymmetric positive definite linear systems. Unlike CG, SCG requires the solution of a lower triangular linear system to produce each semi-conjugate direction. We prove that SCG is theoretically equivalent to the full orthogonalization method (FOM), which is based on the Arnoldi process and converges in a finite number of steps. Because SCG's triangular system increases in size each iteration, we study a sliding window implementation (SWI) to improve efficiency, and show that the directions produced are still locally semi-conjugate. A counterexample illustrates that SWI is different from the direct incomplete orthogonalization method (DIOM), which is FOM with a sliding window. Numerical experiments from the convection-diffusion equation and other applications show that SCG is robust and that the sliding window implementation SWI allows SCG to solve large systems efficiently.

We consider the generalized successive overrelaxation (GSOR) method for solving a class of block three-by-three saddle-point problems. Based on the necessary and sufficient conditions for all roots of a real cubic polynomial to have modulus less than one, we derive convergence results under reasonable assumptions. We also analyze a class of block lower triangular preconditioners induced from GSOR and derive explicit and sharp spectral bounds for the preconditioned matrices. We report numerical experiments on test problems from the liquid crystal director model and the coupled Stokes-Darcy flow, demonstrating the usefulness of GSOR.

We propose iterative projection methods for solving square or rectangular consistent linear systems Ax = b. Existing projection methods use sketching matrices (possibly randomized) to generate a sequence of small projected subproblems, but even the smaller systems can be costly. We develop a process that appends one column to the sketching matrix each iteration and converges in a finite number of iterations whether the sketch is random or deterministic. In general, our process generates orthogonal updates to the approximate solution xk. By choosing the sketch to be the set of all previous residuals, we obtain a simple recursive update and convergence in at most rank(A) iterations (in exact arithmetic). By choosing a sequence of identity columns for the sketch, we develop a generalization of the Kaczmarz method. In experiments on large sparse systems, our method (PLSS) with residual sketches is competitive with LSQR and LSMR, and with residual and identity sketches compares favorably with state-of-the-art randomized methods.

We propose a solution strategy for linear systems arising in interior method optimization, which is suitable for implementation on hardware accelerators such as graphical processing units (GPUs). The current gold standard for solving these systems is the LDL^T factorization. However, LDL^T requires pivoting during factorization, which substantially increases communication cost and degrades performance on GPUs. Our novel approach solves a large indefinite system by solving multiple smaller positive definite systems, using an iterative solve for the Schur complement and an inner direct solve (via Cholesky factorization) within each iteration. Cholesky is stable without pivoting, thereby reducing communication and allowing reuse of the symbolic factorization. We demonstrate the practicality of our approach and show that on large systems it can efficiently utilize GPUs and outperform LDL^T factorization of the full system.

Abstract We delve into the fascinating crossroads of artificial intelligence (AI) and cognitive science, spotlighting the OpenAI advanced language model, ChatGPT. Renowned for generating human-like text, ChatGPT has been widely used in various applications. However, its ability to replicate human cognitive processes, particularly decision-making behavior, remains largely unexplored and untapped. We evaluate ChatGPT's decision-making patterns and show that they strikingly mirror those of human subjects, even patterns traditionally termed ''irrational'' under standard economic theory. This finding challenges the prevailing assumption that AI systems operate solely on rational computations. It suggests that, despite its algorithmic nature, ChatGPT can reflect human cognitive biases when simulating human roles, thus adding a new dimension to our understanding of AI behaviour. Our result places AI models like ChatGPT in a broader context of cognitive science, indicating their potential to mimic not just human language but also human cognitive processes. From a broader perspective, our findings underscore the capacity of AI in behavioral research and stimulate a necessary dialogue on AI design, transparency, and ethical implications. Our study bridges human and machine intelligence, highlighting the potential of AI to enhance our understanding of decision-making processes in artificial agents.

We provide eigenvalues bounds for a new formulation of the step equations in interior methods for convex quadratic optimization. The matrix of our formulation, named K-2.5, has bounded condition number, converges to a well-defined limit under strict complementarity, and has the same size as the traditional, ill-conditioned, saddle-point formulation. We evaluate the performance in the context of a Matlab object-oriented implementation of PDCO, an interior-point solver for minimizing a smooth convex function subject to linear constraints. The main benefit of our implementation, named PDCOO, is to separate the logic of the interior-point method from the formulation of the system used to compute a step at each iteration and the method used to solve the system. Thus, PDCOO allows easy addition of a new system formulation and/or solution method for experimentation. Our numerical experiments indicate that the K-2.5 formulation has the same storage requirements as the traditional ill-conditioned saddle-point formulation, and its condition is often more favourable than the unsymmetric block 3 x 3 formulation.

The linear equations that arise in interior methods for constrained optimization are sparse symmetric indefinite and become extremely ill-conditioned as the interior method converges. These linear systems present a challenge for existing solver frameworks based on sparse LU or LDL^T decompositions. We benchmark five well known direct linear solver packages using matrices extracted from power grid optimization problems. The achieved solution accuracy varies greatly among the packages. None of the tested packages delivers significant GPU acceleration for our test cases.

For optimization problems with linear equality constraints, we prove that the (1,1) block of the inverse KKT matrix remains unchanged when projected onto the nullspace of the constraint matrix. We develop reduced compact representations of the limited-memory inverse BFGS Hessian to compute search directions efficiently when the constraint Jacobian is sparse. Orthogonal projections are implemented by a sparse QR factorization or a preconditioned LSQR iteration. In numerical experiments two proposed trust-region algorithms improve in computation times, often significantly, compared to previous implementations of related algorithms and compared to IPOPT.

Algorithm NCL is designed for general smooth optimization problems where first and second derivatives are available, including problems whose constraints may not be linearly independent at a solution (i.e., do not satisfy the LICQ). It is equivalent to the LANCELOT augmented Lagrangian method, reformulated as a short sequence of nonlinearly constrained subproblems that can be solved efficiently by IPOPT and KNITRO, with warm starts on each subproblem. We give numerical results from a Julia implementation of Algorithm NCL on tax policy models that do not satisfy the LICQ, and on nonlinear least-squares problems and general problems from the CUTEst test set.

Starting with a comprehensive generic reconstruction of human metabolism, we generated high-quality, constraint-based, genome-scale, cell-type and condition specific models of metabolism in human dopaminergic neurons, the cell type most vulnerable to degeneration in Parkinson ’ s disease. They are a synthesis of extensive manual curation of the biochemical literature on neuronal metabolism, together with novel, quantitative, transcriptomic and targeted exometabolomic data from human stem cell-derived, midbrainspecific, dopaminergic neurons in vitro . Thermodynamic constraint-based modelling enabled qualitatively accurate and moderately quantitatively accurate prediction of dopaminergic neuronal metabolite exchange fluxes, including predicting the consequences of metabolic perturbations in a manner also consistent with literature on monogenic mitochondrial Parkinson ’ s disease. These dopaminergic neurons models provide a foundation for a quantitative systems biochemistry approach to metabolic dysfunction in Parkinson ’ s disease. Moreover, the plethora of novel mathematical and computational approaches required to develop them are generalisable to study any other disease associated with metabolic dysfunction.

In time series analysis, when fitting an autoregressive model, one must solve a Toeplitz ordinary least squares problem numerous times to find an appropriate model, which can severely affect computational times with large data sets. Two recent algorithms (LSAR and Repeated Halving) have applied randomized numerical linear algebra (RandNLA) techniques to fitting an autoregressive model to big time-series data. We investigate and compare the quality of these two approximation algorithms on large-scale synthetic and real-world data. While both algorithms display comparable results for synthetic datasets, the LSAR algorithm appears to be more robust when applied to real-world time series data. We conclude that RandNLA is effective in the context of big-data time series.

We recently showed, in a simulation study using two artificial signals, that our PDCO (Primal Dual interior method for Convex Objectives) reconstruction algorithm can be efficiently used for the reconstruction of low-field proton nuclear magnetic resonance (1H LF-NMR) relaxation signals into T1 (spin–lattice) vs. T2 (spin–spin) time 2D graphs of a material’s composition. In the present study, for highly complex materials, we demonstrate the PDCO’s reconstruction efficacy for a much wider range of simulated signals with higher complexity and different signal-to-noise ratios (SNR) taken from actual reconstructed 1H LF-NMR spectroscopy signals of oleic acid and cattle manure. The optimal regularization parameters of the PDCO’s reconstructing algorithm were identified for this large range of simulated LF-NMR signals and noise values. These simulated compact graphical and numerical representations demonstrated 1H LF-NMR relaxation signals of complex materials can be accurately reconstructed into T1 − T2 time graphs of a material’s chemical and morphology. The present study further confirmed that an optimal single set of regulatory parameters for the data reconstruction algorithms could be used for different materials or different batches of the same material.

We develop a general equality-constrained nonlinear optimization algorithm based on a smooth penalty function proposed by Fletcher (1970). Although it was historically considered to be computationally prohibitive in practice, we demonstrate that the computational kernels required are no more expensive than other widely accepted methods for nonlinear optimization. The main kernel required to evaluate the penalty function and its derivatives is solving a structured linear system. We show how to solve this system efficiently by storing a single factorization each iteration when the matrices are available explicitly. We further show how to adapt the penalty function to the class of factorization-free algorithms by solving the linear system iteratively. The penalty function therefore has promise when the linear system can be solved efficiently, e.g., for PDE-constrained optimization problems where efficient preconditioners exist. We discuss extensions including handling simple constraints explicitly, regularizing the penalty function, and inexact evaluation of the penalty function and its gradients. We demonstrate the merits of the approach and its various features on some nonlinear programs from a standard test set, and some PDE-constrained optimization problems.