Michael Saunders
Professor (Research) of Management Science and Engineering, Emeritus
Bio
Saunders develops mathematical methods for solving largescale constrained optimization problems and large systems of equations. He also implements such methods as generalpurpose software to allow their use in many areas of engineering, science, and business. He is codeveloper of the largescale optimizers MINOS, SNOPT, SQOPT, PDCO, the dense QP and NLP solvers LSSOL, QPOPT, NPSOL, and the linear equation solvers SYMMLQ, MINRES, MINRESQLP, LSQR, LSMR, LSLQ, LNLQ, LSRN, LUSOL.
Honors & Awards

OrchardHays Prize, MPS (1985)

Highly Cited Researcher, Computer Science, ISI (2004)

Highly Cited Researcher, Mathematics, ISI (2007)

Honorary Fellow, RSNZ (2007)

Linear Algebra Prize, SIAM (2012)

Invention Hall of Fame, OTL, Stanford University (2012)

Fellow, SIAM (2013)
Boards, Advisory Committees, Professional Organizations

Associate Editor, NACO (2010  2016)

Member, ACM (1982  Present)

Member, INFORMS (2010  Present)

Member, ORSNZ (1990  Present)

Member, SIAM (1980  Present)

Associate Editor, ACM TOMS (1982  2004)

Associate Editor, SIAM Journal on Optimization (1989  2002)

Associate Editor, OPTE (1999  Present)
Professional Education

B.Sc. (Hons), Canterbury, Mathematics (1965)

MS, Stanford University, Computer Science (1970)

PhD, Stanford University, Computer Science (1972)
Stanford Advisees

Doctoral Dissertation Reader (AC)
Laura Lyman 
Doctoral Dissertation CoAdvisor (NonAC)
Shaked Regev
All Publications

IMPLEMENTING A SMOOTH EXACT PENALTY FUNCTION FOR EQUALITYCONSTRAINED NONLINEAR OPTIMIZATION
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2020; 42 (3): A1809–A1835
View details for DOI 10.1137/19M1238265
View details for Web of Science ID 000551255700006

IMPLEMENTING A SMOOTH EXACT PENALTY FUNCTION FOR GENERAL CONSTRAINED NONLINEAR OPTIMIZATION
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2020; 42 (3): A1836–A1859
View details for DOI 10.1137/19M1255069
View details for Web of Science ID 000551255700011

Analysis of the Regularization Parameters of PrimalDual Interior Method for Convex Objectives Applied to H1 Low Field Nuclear Magnetic Resonance Data Processing (vol 49, pg 1129, 2018)
APPLIED MAGNETIC RESONANCE
2019; 50 (13): 521
View details for DOI 10.1007/s0072301810724
View details for Web of Science ID 000458124800037

EUCLIDEANNORM ERROR BOUNDS FOR SYMMLQ AND CG
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
2019; 40 (1): 235–53
View details for DOI 10.1137/16M1094816
View details for Web of Science ID 000462583900011

LSLQ: AN ITERATIVE METHOD FOR LINEAR LEASTSQUARES WITH AN ERROR MINIMIZATION PROPERTY
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
2019; 40 (1): 254–75
View details for DOI 10.1137/17M1113552
View details for Web of Science ID 000462583900012

Stabilized Optimization Via an NCL Algorithm
SPRINGER. 2018: 173–91
View details for DOI 10.1007/9783319900261_8
View details for Web of Science ID 000448498900008

Reliable and efficient solution of genomescale models of Metabolism and macromolecular Expression
SCIENTIFIC REPORTS
2017; 7
Abstract
ConstraintBased Reconstruction and Analysis (COBRA) is currently the only methodology that permits integrated modeling of Metabolism and macromolecular Expression (ME) at genomescale. Linear optimization computes steadystate flux solutions to ME models, but flux values are spread over many orders of magnitude. Data values also have greatly varying magnitudes. Standard doubleprecision solvers may return inaccurate solutions or report that no solution exists. Exact simplex solvers based on rational arithmetic require a nearoptimal warm start to be practical on large problems (current ME models have 70,000 constraints and variables and will grow larger). We have developed a quadrupleprecision version of our linear and nonlinear optimizer MINOS, and a solution procedure (DQQ) involving Double and Quad MINOS that achieves reliability and efficiency for ME models and other challenging problems tested here. DQQ will enable extensive use of large linear and nonlinear models in systems biology and other applications involving multiscale data.
View details for DOI 10.1038/srep40863
View details for Web of Science ID 000392188100001
View details for PubMedID 28098205
View details for PubMedCentralID PMC5241643

Conditions for duality between fluxes and concentrations in biochemical networks
JOURNAL OF THEORETICAL BIOLOGY
2016; 409: 110
Abstract
Mathematical and computational modelling of biochemical networks is often done in terms of either the concentrations of molecular species or the fluxes of biochemical reactions. When is mathematical modelling from either perspective equivalent to the other? Mathematical duality translates concepts, theorems or mathematical structures into other concepts, theorems or structures, in a onetoone manner. We present a novel stoichiometric condition that is necessary and sufficient for duality between unidirectional fluxes and concentrations. Our numerical experiments, with computational models derived from a range of genomescale biochemical networks, suggest that this fluxconcentration duality is a pervasive property of biochemical networks. We also provide a combinatorial characterisation that is sufficient to ensure fluxconcentration duality.The condition prescribes that, for every two disjoint sets of molecular species, there is at least one reaction complex that involves species from only one of the two sets. When unidirectional fluxes and molecular species concentrations are dual vectors, this implies that the behaviour of the corresponding biochemical network can be described entirely in terms of either concentrations or unidirectional fluxes.
View details for DOI 10.1016/j.jtbi.2016.06.033
View details for Web of Science ID 000385471800001
View details for PubMedID 27345817
View details for PubMedCentralID PMC5048525
 Laplace inversion of lowresolution NMR relaxometry data using sparse representation methods Concepts in Magnetic Resonance Part A 2013; 42A:3: 7288
 Novel 1H low field nuclear magnetic resonance applications for the field of biodiesel Biotechnologyfor Biofuels 2013; 6:55: 20
 LSRN: a parallel iterative solver for strongly over or underdetermined systems SIAM J. Sci. Comp. 2013; 36 (2): C95C118

A variational principle for computing nonequilibrium fluxes and potentials in genomescale biochemical networks
JOURNAL OF THEORETICAL BIOLOGY
2012; 292: 7177
Abstract
We derive a convex optimization problem on a steadystate nonequilibrium network of biochemical reactions, with the property that energy conservation and the second law of thermodynamics both hold at the problem solution. This suggests a new variational principle for biochemical networks that can be implemented in a computationally tractable manner. We derive the Lagrange dual of the optimization problem and use strong duality to demonstrate that a biochemical analogue of Tellegen's theorem holds at optimality. Each optimal flux is dependent on a free parameter that we relate to an elementary kinetic parameter when mass action kinetics is assumed.
View details for DOI 10.1016/j.jtbi.2011.09.029
View details for Web of Science ID 000297450100008
View details for PubMedID 21983269

LSMR: AN ITERATIVE ALGORITHM FOR SPARSE LEASTSQUARES PROBLEMS
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2011; 33 (5): 29502971
View details for Web of Science ID 000296591200039
 SNOPT: An SQP algorithm for largescaleconstrained optimization, SIGEST article SIAM Rev. 2005; 1 (47): 99131
 Atomic decomposition by basis pursuit, SIGEST article SIAM Rev. 2001; 1 (43): 129159

HyKKT: a hybrid directiterative method for solving KKT linear systems
OPTIMIZATION METHODS & SOFTWARE
2022
View details for DOI 10.1080/10556788.2022.2124990
View details for Web of Science ID 000860640100001

Linear solvers for power grid optimization problems: A review of GPUaccelerated linear solvers
PARALLEL COMPUTING
2022; 111
View details for DOI 10.1016/j.parco.2021.102870
View details for Web of Science ID 000800000400004

LARGESCALE OPTIMIZATION WITH LINEAR EQUALITY CONSTRAINTS USING REDUCED COMPACT REPRESENTATION\ast
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2022; 44 (1): A103A127
View details for DOI 10.1137/21M1393819
View details for Web of Science ID 000773632300005

Linear systems arising in interior methods for convex optimization: a symmetric formulation with bounded condition number
OPTIMIZATION METHODS & SOFTWARE
2021
View details for DOI 10.1080/10556788.2021.1965599
View details for Web of Science ID 000705433600001

SimulationBased Sensitivity Analysis of Regularization Parameters for Robust Reconstruction of Complex Material's T1  (T2H)H1 LFNMR Energy Relaxation Signals
APPLIED MAGNETIC RESONANCE
2019
View details for DOI 10.1007/s00723019011731
View details for Web of Science ID 000497186200001

Creation and analysis of biochemical constraintbased models using the COBRA Toolbox v.3.0.
Nature protocols
2019
Abstract
Constraintbased reconstruction and analysis (COBRA) provides a molecular mechanistic framework for integrative analysis of experimental molecular systems biology data and quantitative prediction of physicochemically and biochemically feasible phenotypic states. The COBRA Toolbox is a comprehensive desktop software suite of interoperable COBRA methods. It has found widespread application in biology, biomedicine, and biotechnology because its functions can be flexibly combined to implement tailored COBRA protocols for any biochemical network. This protocol is an update to the COBRA Toolbox v.1.0 and v.2.0. Version 3.0 includes new methods for qualitycontrolled reconstruction, modeling, topological analysis, strain and experimental design, and network visualization, as well as network integration of chemoinformatic, metabolomic, transcriptomic, proteomic, and thermochemical data. New multilingual code integration also enables an expansion in COBRA application scope via highprecision, highperformance, and nonlinear numerical optimization solvers for multiscale, multicellular, and reaction kinetic modeling, respectively. This protocol provides an overview of all these new features and can be adapted to generate and analyze constraintbased models in a wide variety of scenarios. The COBRA Toolbox v.3.0 provides an unparalleled depth of COBRA methods.
View details for PubMedID 30787451

DynamicME: dynamic simulation and refinement of integrated models of metabolism and protein expression.
BMC systems biology
2019; 13 (1): 2
Abstract
BACKGROUND: Genomescale models of metabolism and macromolecular expression (ME models) enable systemslevel computation of proteome allocation coupled to metabolic phenotype.RESULTS: We develop DynamicME, an algorithm enabling timecourse simulation of cell metabolism and protein expression. DynamicME correctly predicted the substrate utilization hierarchy on a mixed carbon substrate medium. We also found good agreement between predicted and measured timecourse expression profiles. ME models involve considerably more parameters than metabolic models (M models). We thus generate an ensemble of models (each model having its rate constants perturbed), and then analyze the models by identifying archetypal timecourse metabolite concentration profiles. Furthermore, we use a metaheuristic optimization method to calibrate ME model parameters using timecourse measurements such as from a (fed) batch culture. Finally, we show that constraints on protein concentration dynamics ("inertia") alter the metabolic response to environmental fluctuations, including increased substratelevel phosphorylation and lowered oxidative phosphorylation.CONCLUSIONS: Overall, DynamicME provides a novel method for understanding proteome allocation and metabolism under complex and transient environments, and to utilize timecourse cell culture data for modelbased interpretation or model refinement.
View details for PubMedID 30626386

LNLQ: AN ITERATIVE METHOD FOR LEASTNORM PROBLEMS WITH AN ERROR MINIMIZATION PROPERTY
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
2019; 40 (3): 1102–24
View details for DOI 10.1137/18M1194948
View details for Web of Science ID 000487856100012

Estimating Cellular Goals from HighDimensional Biological Data
ASSOC COMPUTING MACHINERY. 2019: 2202–11
View details for DOI 10.1145/3292500.3330775
View details for Web of Science ID 000485562502026

Analysis of the Regularization Parameters of PrimalDual Interior Method for Convex Objectives Applied to H1 Low Field Nuclear Magnetic Resonance Data Processing
APPLIED MAGNETIC RESONANCE
2018; 49 (10): 1129–50
View details for DOI 10.1007/s0072301810484
View details for Web of Science ID 000444236900007

Principles of proteome allocation are revealed using proteomic data and genomescale models
SCIENTIFIC REPORTS
2016; 6
Abstract
Integrating omics data to refine or make contextspecific models is an active field of constraintbased modeling. Proteomics now cover over 95% of the Escherichia coli proteome by mass. Genomescale models of Metabolism and macromolecular Expression (ME) compute proteome allocation linked to metabolism and fitness. Using proteomics data, we formulated allocation constraints for key proteome sectors in the ME model. The resulting calibrated model effectively computed the "generalist" (wildtype) E. coli proteome and phenotype across diverse growth environments. Across 15 growth conditions, prediction errors for growth rate and metabolic fluxes were 69% and 14% lower, respectively. The sectorconstrained ME model thus represents a generalist ME model reflecting both growth rate maximization and "hedging" against uncertain environments and stresses, as indicated by significant enrichment of these sectors for the general stress response sigma factor σ(S). Finally, the sector constraints represent a general formalism for integrating omics data from any experimental condition into constraintbased ME models. The constraints can be finegrained (individual proteins) or coarsegrained (functionallyrelated protein groups) as demonstrated here. This flexible formalism provides an accessible approach for narrowing the gap between the complexity captured by omics data and governing principles of proteome allocation described by systemslevel models.
View details for DOI 10.1038/srep36734
View details for Web of Science ID 000388074400001
View details for PubMedID 27857205
View details for PubMedCentralID PMC5114563

solveME: fast and reliable solution of nonlinear ME models.
BMC bioinformatics
2016; 17 (1): 391
Abstract
Genomescale models of metabolism and macromolecular expression (ME) significantly expand the scope and predictive capabilities of constraintbased modeling. ME models present considerable computational challenges: they are much (>30 times) larger than corresponding metabolic reconstructions (M models), are multiscale, and growth maximization is a nonlinear programming (NLP) problem, mainly due to macromolecule dilution constraints.Here, we address these computational challenges. We develop a fast and numerically reliable solution method for growth maximization in ME models using a quadprecision NLP solver (Quad MINOS). Our method was up to 45 % faster than binary search for six significant digits in growth rate. We also develop a fast, quadprecision flux variability analysis that is accelerated (up to 60× speedup) via solver warmstarts. Finally, we employ the tools developed to investigate growthcoupled succinate overproduction, accounting for proteome constraints.Just as genomescale metabolic reconstructions have become an invaluable tool for computational and systems biologists, we anticipate that these fast and numerically reliable ME solution methods will accelerate the widespread adoption of ME models for researchers in these fields.
View details for DOI 10.1186/s1285901612401
View details for PubMedID 27659412
View details for PubMedCentralID PMC5034503

Heart rate analysis by sparse representation for acute pain detection
MEDICAL & BIOLOGICAL ENGINEERING & COMPUTING
2016; 54 (4): 595606
Abstract
Objective pain assessment methods pose an advantage over the currently used subjective pain rating tools. Advanced signal processing methodologies, including the wavelet transform (WT) and the orthogonal matching pursuit algorithm (OMP), were developed in the past two decades. The aim of this study was to apply and compare these timespecific methods to heart rate samples of healthy subjects for acute pain detection. Fifteen adult volunteers participated in a study conducted in the pain clinic at a single center. Each subject's heart rate was sampled for 5min baseline, followed by a cold pressor test (CPT). Analysis was done by the WT and the OMP algorithm with a Fourier/Wavelet dictionary separately. Data from 11 subjects were analyzed. Compared to baseline, The WT analysis showed a significant coefficients' density increase during the pain incline period (p < 0.01) and the entire CPT (p < 0.01), with significantly higher coefficient amplitudes. The OMP analysis showed a significant wavelet coefficients' density increase during pain incline and decline periods (p < 0.01, p < 0.05) and the entire CPT (p < 0.001), with suggestive higher amplitudes. Comparison of both methods showed that during the baseline there was a significant reduction in wavelet coefficient density using the OMP algorithm (p < 0.001). Analysis by the twoway ANOVA with repeated measures showed a significant proportional increase in wavelet coefficients during the incline period and the entire CPT using the OMP algorithm (p < 0.01). Both methods provided accurate and nondelayed detection of pain events. Statistical analysis proved the OMP to be by far more specific allowing the Fourier coefficients to represent the signal's basic harmonics and the wavelet coefficients to focus on the timespecific painful event. This is an initial study using OMP for pain detection; further studies need to prove the efficiency of this system in different settings.
View details for DOI 10.1007/s1151701513503
View details for Web of Science ID 000373021100004

Heart rate analysis by sparse representation for acute pain detection.
Medical & biological engineering & computing
2016; 54 (4): 595606
Abstract
Objective pain assessment methods pose an advantage over the currently used subjective pain rating tools. Advanced signal processing methodologies, including the wavelet transform (WT) and the orthogonal matching pursuit algorithm (OMP), were developed in the past two decades. The aim of this study was to apply and compare these timespecific methods to heart rate samples of healthy subjects for acute pain detection. Fifteen adult volunteers participated in a study conducted in the pain clinic at a single center. Each subject's heart rate was sampled for 5min baseline, followed by a cold pressor test (CPT). Analysis was done by the WT and the OMP algorithm with a Fourier/Wavelet dictionary separately. Data from 11 subjects were analyzed. Compared to baseline, The WT analysis showed a significant coefficients' density increase during the pain incline period (p < 0.01) and the entire CPT (p < 0.01), with significantly higher coefficient amplitudes. The OMP analysis showed a significant wavelet coefficients' density increase during pain incline and decline periods (p < 0.01, p < 0.05) and the entire CPT (p < 0.001), with suggestive higher amplitudes. Comparison of both methods showed that during the baseline there was a significant reduction in wavelet coefficient density using the OMP algorithm (p < 0.001). Analysis by the twoway ANOVA with repeated measures showed a significant proportional increase in wavelet coefficients during the incline period and the entire CPT using the OMP algorithm (p < 0.01). Both methods provided accurate and nondelayed detection of pain events. Statistical analysis proved the OMP to be by far more specific allowing the Fourier coefficients to represent the signal's basic harmonics and the wavelet coefficients to focus on the timespecific painful event. This is an initial study using OMP for pain detection; further studies need to prove the efficiency of this system in different settings.
View details for DOI 10.1007/s1151701513503
View details for PubMedID 26264057

A Practical Factorization of a Schur Complement for PDEConstrained Distributed Optimal Control
JOURNAL OF SCIENTIFIC COMPUTING
2015; 65 (2): 576597
View details for DOI 10.1007/s1091501499760
View details for Web of Science ID 000362911900007

Do genomescale models need exact solvers or clearer standards?
MOLECULAR SYSTEMS BIOLOGY
2015; 11 (10): 831
View details for PubMedID 26467284

Systems biology definition of the core proteome of metabolism and expression is consistent with highthroughput data.
Proceedings of the National Academy of Sciences of the United States of America
2015; 112 (34): 1081010815
Abstract
Finding the minimal set of gene functions needed to sustain life is of both fundamental and practical importance. Minimal gene lists have been proposed by using comparative genomicsbased core proteome definitions. A definition of a core proteome that is supported by empirical data, is understood at the systemslevel, and provides a basis for computing essential cell functions is lacking. Here, we use a systems biologybased genomescale model of metabolism and expression to define a functional core proteome consisting of 356 gene products, accounting for 44% of the Escherichia coli proteome by mass based on proteomics data. This systems biology core proteome includes 212 genes not found in previous comparative genomicsbased core proteome definitions, accounts for 65% of known essential genes in E. coli, and has 78% gene function overlap with minimal genomes (Buchnera aphidicola and Mycoplasma genitalium). Based on transcriptomics data across environmental and genetic backgrounds, the systems biology core proteome is significantly enriched in nondifferentially expressed genes and depleted in differentially expressed genes. Compared with the noncore, core gene expression levels are also similar across genetic backgrounds (two times higher Spearman rank correlation) and exhibit significantly more complex transcriptional and posttranscriptional regulatory features (40% more transcription start sites per gene, 22% longer 5'UTR). Thus, genomescale systems biology approaches rigorously identify a functional core proteome needed to support growth. This framework, validated by using highthroughput datasets, facilitates a mechanistic understanding of systemslevel core proteome function through in silico models; it de facto defines a paleome.
View details for DOI 10.1073/pnas.1501384112
View details for PubMedID 26261351

Study of liquidphase molecular packing interactions and morphology of fatty acid methyl esters (biodiesel).
Biotechnology for biofuels
2015; 8: 12?
Abstract
(1)H low field nuclear magnetic resonance (LFNMR) relaxometry has been suggested as a tool to distinguish between different molecular ensembles in complex systems with differential segmental or whole molecular motion and/or different morphologies. In biodiesel applications the molecular structure versus liquidphase packing morphologies of fatty acid methyl esters (FAMEs) influences physicochemical characteristics of the fuel, including flow properties, operability during cold weather, blending, and more. Still, their liquid morphological structures have scarcely been studied. It was therefore the objective of this work to explore the potential of this technology for characterizing the molecular organization of FAMEs in the liquid phase. This was accomplished by using a combination of supporting advanced technologies.We show that pure oleic acid (OA) and methyl oleate (MO) standards exhibited both similarities and differences in the (1)H LFNMR relaxation times (T2s) and peak areas, for a range of temperatures. Based on Xray measurements, both molecules were found to possess a liquid crystallike order, although a larger fluidity was found for MO, because as the temperature is increased, MO molecules separate both longitudinally and transversely from one another. In addition, both molecules exhibited a preferred direction of diffusion based on the apparent hydrodynamic radius. The close molecular packing arrangement and interactions were found to affect the translational and segmental motions of the molecules, as a result of dimerization of the head group in OA as opposed to weaker polar interactions in MO.A comprehensive model for the liquid crystallike arrangement of FAMEs in the liquid phase is suggested. The differences in translational and segmental motions of the molecules were rationalized by the differences in the (1)H LFNMR T2 distributions of OA and MO, which was further supported by (13)C high field (HF)NMR spectra and (1)H HFNMR relaxation. The proposed assignment allows for material characterization based on parameters that contribute to properties in applications such as biodiesel fuels.
View details for DOI 10.1186/s1306801401947
View details for PubMedID 25688289
View details for PubMedCentralID PMC4329664

LSRN: A PARALLEL ITERATIVE SOLVER FOR STRONGLY OVER OR UNDERDETERMINED SYSTEMS.
SIAM journal on scientific computing : a publication of the Society for Industrial and Applied Mathematics
2014; 36 (2): C95C118
Abstract
We describe a parallel iterative least squares solver named LSRN that is based on random normal projection. LSRN computes the minlength solution to min x∈ℝ n ‖Ax  b‖2, where A ∈ ℝ m × n with m ≫ n or m ≪ n, and where A may be rankdeficient. Tikhonov regularization may also be included. Since A is involved only in matrixmatrix and matrixvector multiplications, it can be a dense or sparse matrix or a linear operator, and LSRN automatically speeds up when A is sparse or a fast linear operator. The preconditioning phase consists of a random normal projection, which is embarrassingly parallel, and a singular value decomposition of size ⌈γ min(m, n)⌉ × min(m, n), where γ is moderately larger than 1, e.g., γ = 2. We prove that the preconditioned system is wellconditioned, with a strong concentration result on the extreme singular values, and hence that the number of iterations is fully predictable when we apply LSQR or the Chebyshev semiiterative method. As we demonstrate, the Chebyshev method is particularly efficient for solving large problems on clusters with high communication cost. Numerical results show that on a sharedmemory machine, LSRN is very competitive with LAPACK's DGELSD and a fast randomized least squares solver called Blendenpik on large dense problems, and it outperforms the least squares solver from SuiteSparseQR on sparse problems without sparsity patterns that can be exploited to reduce fillin. Further experiments show that LSRN scales well on an Amazon Elastic Compute Cloud cluster.
View details for DOI 10.1137/120866580
View details for PubMedID 25419094
View details for PubMedCentralID PMC4238893

Algorithm 937: MINRESQLP for Symmetric and Hermitian Linear Equations and LeastSquares Problems
ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
2014; 40 (2)
View details for DOI 10.1145/2527267
View details for Web of Science ID 000333653400008

Algorithm 937: MINRESQLP for Symmetric and Hermitian Linear Equations and LeastSquares Problems.
ACM transactions on mathematical software. Association for Computing Machinery
2014; 40 (2)
Abstract
We describe algorithm MINRESQLP and its FORTRAN 90 implementation for solving symmetric or Hermitian linear systems or leastsquares problems. If the system is singular, MINRESQLP computes the unique minimumlength solution (also known as the pseudoinverse solution), which generally eludes MINRES. In all cases, it overcomes a potential instability in the original MINRES algorithm. A positivedefinite preconditioner may be supplied. Our FORTRAN 90 implementation illustrates a design pattern that allows users to make problem data known to the solver but hidden and secure from other program units. In particular, we circumvent the need for reverse communication. Example test programs input and solve real or complex problems specified in Matrix Market format. While we focus here on a FORTRAN 90 implementation, we also provide and maintain MATLAB versions of MINRES and MINRESQLP.
View details for DOI 10.1145/2527267
View details for PubMedID 25328255
View details for PubMedCentralID PMC4199394

LSRN: A PARALLEL ITERATIVE SOLVER FOR STRONGLY OVER OR UNDERDETERMINED SYSTEMS
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2014; 36 (2): C95C118
Abstract
We describe a parallel iterative least squares solver named LSRN that is based on random normal projection. LSRN computes the minlength solution to min x∈ℝ n ‖Ax  b‖2, where A ∈ ℝ m × n with m ≫ n or m ≪ n, and where A may be rankdeficient. Tikhonov regularization may also be included. Since A is involved only in matrixmatrix and matrixvector multiplications, it can be a dense or sparse matrix or a linear operator, and LSRN automatically speeds up when A is sparse or a fast linear operator. The preconditioning phase consists of a random normal projection, which is embarrassingly parallel, and a singular value decomposition of size ⌈γ min(m, n)⌉ × min(m, n), where γ is moderately larger than 1, e.g., γ = 2. We prove that the preconditioned system is wellconditioned, with a strong concentration result on the extreme singular values, and hence that the number of iterations is fully predictable when we apply LSQR or the Chebyshev semiiterative method. As we demonstrate, the Chebyshev method is particularly efficient for solving large problems on clusters with high communication cost. Numerical results show that on a sharedmemory machine, LSRN is very competitive with LAPACK's DGELSD and a fast randomized least squares solver called Blendenpik on large dense problems, and it outperforms the least squares solver from SuiteSparseQR on sparse problems without sparsity patterns that can be exploited to reduce fillin. Further experiments show that LSRN scales well on an Amazon Elastic Compute Cloud cluster.
View details for DOI 10.1137/120866580
View details for Web of Science ID 000335817600030
View details for PubMedCentralID PMC4238893

PROXIMAL NEWTONTYPE METHODS FOR MINIMIZING COMPOSITE FUNCTIONS
SIAM JOURNAL ON OPTIMIZATION
2014; 24 (3): 14201443
View details for DOI 10.1137/130921428
View details for Web of Science ID 000343229000019

Robust flux balance analysis of multiscale biochemical reaction networks
BMC BIOINFORMATICS
2013; 14
Abstract
Biological processes such as metabolism, signaling, and macromolecular synthesis can be modeled as large networks of biochemical reactions. Large and comprehensive networks, like integrated networks that represent metabolism and macromolecular synthesis, are inherently multiscale because reaction rates can vary over many orders of magnitude. They require special methods for accurate analysis because naive use of standard optimization systems can produce inaccurate or erroneously infeasible results.We describe techniques enabling offtheshelf optimization software to compute accurate solutions to the poorly scaled optimization problems arising from flux balance analysis of multiscale biochemical reaction networks. We implement lifting techniques for flux balance analysis within the openCOBRA toolbox and demonstrate our techniques using the first integrated reconstruction of metabolism and macromolecular synthesis for E. coli.Our techniques enable accurate flux balance analysis of multiscale networks using offtheshelf optimization software. Although we describe lifting techniques in the context of flux balance analysis, our methods can be used to handle a variety of optimization problems arising from analysis of multiscale network reconstructions.
View details for DOI 10.1186/1471210514240
View details for Web of Science ID 000322915900001
View details for PubMedID 23899245
View details for PubMedCentralID PMC3750362

Equispaced Pareto front construction for constrained biobjective optimization
MATHEMATICAL AND COMPUTER MODELLING
2013; 57 (910): 21222131
View details for DOI 10.1016/j.mcm.2010.12.044
View details for Web of Science ID 000317262100010
 Robust flux balance analysis of multiscale biochemical reaction networks BMC Bioinformatics 2013; 14:240: 6
 CG versus MINRES: An empirical comparison SQUJournal for Science 2012; 17:1: 4462

A HigherOrder Generalized Singular Value Decomposition for Comparison of Global mRNA Expression from Multiple Organisms
PLOS ONE
2011; 6 (12)
Abstract
The number of highdimensional datasets recording multiple aspects of a single phenomenon is increasing in many areas of science, accompanied by a need for mathematical frameworks that can compare multiple largescale matrices with different row dimensions. The only such framework to date, the generalized singular value decomposition (GSVD), is limited to two matrices. We mathematically define a higherorder GSVD (HO GSVD) for N≥2 matrices D(i)∈R(m(i) × n), each with full column rank. Each matrix is exactly factored as D(i)=U(i)Σ(i)V(T), where V, identical in all factorizations, is obtained from the eigensystem SV=VΛ of the arithmetic mean S of all pairwise quotients A(i)A(j)(1) of the matrices A(i)=D(i)(T)D(i), i≠j. We prove that this decomposition extends to higher orders almost all of the mathematical properties of the GSVD. The matrix S is nondefective with V and Λ real. Its eigenvalues satisfy λ(k)≥1. Equality holds if and only if the corresponding eigenvector v(k) is a right basis vector of equal significance in all matrices D(i) and D(j), that is σ(i,k)/σ(j,k)=1 for all i and j, and the corresponding left basis vector u(i,k) is orthogonal to all other vectors in U(i) for all i. The eigenvalues λ(k)=1, therefore, define the "common HO GSVD subspace." We illustrate the HO GSVD with a comparison of genomescale cellcycle mRNA expression from S. pombe, S. cerevisiae and human. Unlike existing algorithms, a mapping among the genes of these disparate organisms is not required. We find that the approximately common HO GSVD subspace represents the cellcycle mRNA expression oscillations, which are similar among the datasets. Simultaneous reconstruction in the common subspace, therefore, removes the experimental artifacts, which are dissimilar, from the datasets. In the simultaneous sequenceindependent classification of the genes of the three organisms in this common subspace, genes of highly conserved sequences but significantly different cellcycle peak times are correctly classified.
View details for DOI 10.1371/journal.pone.0028072
View details for Web of Science ID 000299684700003
View details for PubMedID 22216090

MINRESQLP: A KRYLOV SUBSPACE METHOD FOR INDEFINITE OR SINGULAR SYMMETRIC SYSTEMS
SIAM JOURNAL ON SCIENTIFIC COMPUTING
2011; 33 (4): 18101836
View details for DOI 10.1137/100787921
View details for Web of Science ID 000294293200016

Nonconservative Robust Control: Optimized and Constrained Sensitivity Functions
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY
2009; 17 (2): 298308
View details for DOI 10.1109/TCST.2008.924564
View details for Web of Science ID 000263832000004

STABILIZING POLICY IMPROVEMENT FOR LARGESCALE INFINITEHORIZON DYNAMIC PROGRAMMING
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
2009; 31 (2): 434459
View details for DOI 10.1137/060653305
View details for Web of Science ID 000267745500012

Variational Bayesian image restoration based on a product of tdistributions image prior
IEEE TRANSACTIONS ON IMAGE PROCESSING
2008; 17 (10): 17951805
Abstract
Image priors based on products have been recognized to offer many advantages because they allow simultaneous enforcement of multiple constraints. However, they are inconvenient for Bayesian inference because it is hard to find their normalization constant in closed form. In this paper, a new Bayesian algorithm is proposed for the image restoration problem that bypasses this difficulty. An image prior is defined by imposing Studentt densities on the outputs of local convolutional filters. A variational methodology, with a constrained expectation step, is used to infer the restored image. Numerical experiments are shown that compare this methodology to previous ones and demonstrate its advantages.
View details for DOI 10.1109/TIP.2008.2002828
View details for Web of Science ID 000259372100005
View details for PubMedID 18784028

George B. Dantzig and systems optimization
DISCRETE OPTIMIZATION
2008; 5 (2): 151158
View details for DOI 10.1016/j.disopt.2007.01.002
View details for Web of Science ID 000255475400002

Discussion: The Dantzig selector: Statistical estimation when p is much larger than n
ANNALS OF STATISTICS
2007; 35 (6): 23852391
View details for DOI 10.1214/00905360700000479
View details for Web of Science ID 000253077800007
 Commentary on Methods for modifying matrix factorizations Milestones in Matrix Computation: Selected Works of Gene H. Golub With Commentaries edited by Chan, R., H., Greif, C., O'Leary, D., P. Oxford University Press. 2007: 310–310

SpaseLoc: An adaptive subproblem algorithm for scalable wireless sensor network localization
SIAM JOURNAL ON OPTIMIZATION
2006; 17 (4): 11021128
View details for DOI 10.1137/040621600
View details for Web of Science ID 000244631800007

SNOPT: An SQP algorithm for largescale constrained optimization (Reprinted from SIAM Journal Optimization, vol 12, pg 9791006, 2002)
SIAM REVIEW
2005; 47 (1): 99131
View details for DOI 10.1137/S0036144504446096
View details for Web of Science ID 000227119200005

A globally convergent linearly constrained Lagrangian method for nonlinear optimization
SIAM JOURNAL ON OPTIMIZATION
2005; 15 (3): 863897
View details for DOI 10.1137/S1052623402419789
View details for Web of Science ID 000229826800011

Sparsity and smoothness via the fused lasso
JOURNAL OF THE ROYAL STATISTICAL SOCIETY SERIES BSTATISTICAL METHODOLOGY
2005; 67: 91108
View details for Web of Science ID 000225686900006

A bisection algorithm for the mixed mu upper bound and its supremum
American Control Conference
IEEE. 2004: 2665–2670
View details for Web of Science ID 000224688300453

Subspace preconditioned LSQR for discrete illposed problems
Conference on Computational Linear Algebra with Applications
SPRINGER. 2003: 975–89
View details for Web of Science ID 000188719300011

SNOPT: An SQP algorithm for largescale constrained optimization
SIAM JOURNAL ON OPTIMIZATION
2002; 12 (4): 9791006
View details for Web of Science ID 000175810600007
 Global controller optimization using Horowitz bounds 2002

Atomic decomposition by basis pursuit
SIAM REVIEW
2001; 43 (1): 129159
View details for Web of Science ID 000167366100008

Atomic decomposition by basis pursuit
SIAM JOURNAL ON SCIENTIFIC COMPUTING
1998; 20 (1): 3361
View details for Web of Science ID 000075434800003
 SNOPT: A Fortran software package to solve largescale optimization problems 1998

OSSE mapping of galactic 511 keV positron annihilation line emission
ASTROPHYSICAL JOURNAL
1997; 491 (2): 725748
View details for Web of Science ID 000071152600025

Computing projections with LSQR
BIT NUMERICAL MATHEMATICS
1997; 37 (1): 96104
View details for Web of Science ID A1997WK05500008

Nonparametric estimates of high energy gammaray source distributions
4th Compton Symposium
AIP PRESS. 1997: 1601–5
View details for Web of Science ID 000071400800240

Choleskybased methods for sparse least squares: The benefits of regularization
AMS/IMS/SIAM Summer Research Conference on Linear and Nonlinear Conjugate GradientRelated Methods
SIAM. 1996: 92–100
View details for Web of Science ID A1996BF52D00008
 SQP methods for largescale optimization 1996

On the stability of Cholesky factorization for symmetric quasidefinite systems
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
1996; 17 (1): 3546
View details for Web of Science ID A1996TV10700002

Solution of sparse rectangular systems using LSQR and Craig
BIT NUMERICAL MATHEMATICS
1995; 35 (4): 588604
View details for Web of Science ID A1995TL73600010

Primaldual methods for linear programming
MATHEMATICAL PROGRAMMING
1995; 70 (3): 251277
View details for Web of Science ID A1995TK93600002

A PRACTICAL INTERIORPOINT METHOD FOR CONVEXPROGRAMMING
SIAM JOURNAL ON OPTIMIZATION
1995; 5 (1): 149171
View details for Web of Science ID A1995QJ02000008
 MINOS(IIS) version 4.2: Analyzing infeasibilities inlinear programming Eur. J. Oper. Res. 1995; 81: 217218
 Fortran software for optimization 1995

THE SIMPLEX ALGORITHM WITH A NEW PRIMAL AND DUAL PIVOT RULE
OPERATIONS RESEARCH LETTERS
1994; 16 (3): 121127
View details for Web of Science ID A1994PV30700001

SOLVING REDUCED KKT SYSTEMS IN BARRIER METHODS FOR LINEARPROGRAMMING
15th Dundee Conference on Numerical Analysis
LONGMAN SCIENTIFIC & TECHNICAL. 1994: 89–104
View details for Web of Science ID A1994BA91K00006

Fortran software for optimization
1995 NSF Design and Manufacturing Grantees Conference
SOC MANUFACTURING ENGINEERS. 1994: 31–32
View details for Web of Science ID A1994BG45A00016
 Largescale SQP methods and their applicationin trajectory optimization Control Applications of Optimization edited by Bulirsch, R., Kraft, D. Birkhauser Verlag, Basel,Boston, Stuttgart. 1994: 29–42
 Solving reduced KKT systems in barrier methods for linear programming Numerical Analysis 1993 edited by Watson, G., A., Grffiths, D. Pitman Research Notes in Mathematics 303, Longmans Press. 1994: 89–104
 Major Cholesky would feel proud ORSA J. Comput. 1994; 6: 2327

PRECONDITIONERS FOR INDEFINITE SYSTEMS ARISING IN OPTIMIZATION
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
1992; 13 (1): 292311
View details for Web of Science ID A1992HC84600022
 Some theoretical properties of an augmented Lagrangian merit function Advances in Optimization and Parallel Computing edited by Pardalos, P., M. NorthHolland, Amsterdam. 1992: 101–128
 The applicationof nonlinear programming and collocation to optimal aeroassisted orbital transfers 1992

A BLOCKLU UPDATE FOR LARGESCALE LINEARPROGRAMMING
SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS
1992; 13 (1): 191201
View details for Web of Science ID A1992HC84600016

INERTIACONTROLLING METHODS FOR GENERAL QUADRATICPROGRAMMING
SIAM REVIEW
1991; 33 (1): 136
View details for Web of Science ID A1991FD40400001
 An adaptive primaldual method for linear programming Math.Prog. Soc., Committee on Algorithms Newsletter 1991; 19: 716
 A Schurcomplement method forsparse quadratic programming Reliable Numerical Computation edited by Cox, M., G., Hammarling, S. Oxford University Press, Oxford and New York. 1990: 113–138

A PRACTICAL ANTICYCLING PROCEDURE FOR LINEARLY CONSTRAINED OPTIMIZATION
MATHEMATICAL PROGRAMMING
1989; 45 (3): 437474
View details for Web of Science ID A1989CN43300004
 Constrained nonlinear programming Optimization Handbooks in Operations Research and Management Science edited by Nemhauser, G., L., G., A., H., Kan, R. NorthHolland, Amsterdam. 1989: 171–210

2 CONJUGATEGRADIENTTYPE METHODS FOR UNSYMMETRIC LINEAREQUATIONS
SIAM JOURNAL ON NUMERICAL ANALYSIS
1988; 25 (4): 927940
View details for Web of Science ID A1988P634900009

RECENT DEVELOPMENTS IN CONSTRAINED OPTIMIZATION
JOURNAL OF COMPUTATIONAL AND APPLIED MATHEMATICS
1988; 22 (23): 257270
View details for Web of Science ID A1988P398600009
 Two conjugategradienttype methods forunsymmetric linear equations SIAM J. Numer. Anal. 1988; 25: 927940
 GAMS/MINOS GAMS: A User's Guide edited by Brooke, A., Kendrick, D., Meeraus, A. The Scientic Press. 1988: 201–224

MAINTAINING LU FACTORS OF A GENERAL SPARSEMATRIX
LINEAR ALGEBRA AND ITS APPLICATIONS
1987; 889: 239270
View details for Web of Science ID A1987G801700014

ON PROJECTED NEWTON BARRIER METHODS FOR LINEARPROGRAMMING AND AN EQUIVALENCE TO KARMARKAR PROJECTIVE METHOD
MATHEMATICAL PROGRAMMING
1986; 36 (2): 183209
View details for Web of Science ID A1986F105800006

CONSIDERATIONS OF NUMERICALANALYSIS IN A SEQUENTIAL QUADRATICPROGRAMMING METHOD
LECTURE NOTES IN MATHEMATICS
1986; 1230: 4662
View details for Web of Science ID A1986G659700004
 Considerations of numerical analysis in sequential quadratic programming methods Numerical Analysis edited by Hennart, J., P. SpringerVerlag, New York and London. 1986: 46–62

PROPERTIES OF A REPRESENTATION OF A BASIS FOR THE NULL SPACE
MATHEMATICAL PROGRAMMING
1985; 33 (2): 172186
View details for Web of Science ID A1985ATG1200005
 Software and its relationship tomethods Numerical Optimization 1984 edited by Boggs, P., T., Byrd, R., H., B., R. SIAM, Philadelphia. 1985: 139–159
 Model building and practical aspects of nonlinear programming Computational Mathematical Programming edited by Schittkowski, K. NATO ASI, SpringerVerlag,Berlin and New York. 1985: 209–247

PROCEDURES FOR OPTIMIZATION PROBLEMS WITH A MIXTURE OF BOUNDS AND GENERAL LINEAR CONSTRAINTS
ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
1984; 10 (3): 282298
View details for Web of Science ID A1984TU57100006
 Sequential quadratic programming methods for nonlinear programming Computer Aided Analysis and Optimization of Mechanical System Dynamics edited by Haug, E., J. NATO ASI. 1984: 679–697

SPARSEMATRIX METHODS IN OPTIMIZATION
SIAM JOURNAL ON SCIENTIFIC AND STATISTICAL COMPUTING
1984; 5 (3): 562589
View details for Web of Science ID A1984TG34100006

TRENDS IN NONLINEARPROGRAMMING SOFTWARE
EUROPEAN JOURNAL OF OPERATIONAL RESEARCH
1984; 17 (2): 141149
View details for Web of Science ID A1984TF70000001

AQUIFER RECLAMATION DESIGN  THE USE OF CONTAMINANT TRANSPORT SIMULATION COMBINED WITH NONLINEARPROGRAMMING
WATER RESOURCES RESEARCH
1984; 20 (4): 415427
View details for Web of Science ID A1984SM88600001

A WEIGHTED GRAMSCHMIDT METHOD FOR CONVEX QUADRATICPROGRAMMING
MATHEMATICAL PROGRAMMING
1984; 30 (2): 176195
View details for Web of Science ID A1984TL33800004

COMPUTING FORWARDDIFFERENCE INTERVALS FOR NUMERICAL OPTIMIZATION
SIAM JOURNAL ON SCIENTIFIC AND STATISTICAL COMPUTING
1983; 4 (2): 310321
View details for Web of Science ID A1983QQ77800015

ALGORITHM583  LSQR  SPARSE LINEAREQUATIONS AND LEASTSQUARES PROBLEMS
ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
1982; 8 (2): 195209
View details for Web of Science ID A1982PE14000007
 Software for constrained optimization Nonlinear Optimization 1981 edited by Powell, M. J., D. Academic Press, London and New York. 1982: 381–393
 Linearly constrained optimization Nonlinear Optimization 1981 edited by Powell, M. J., D. Academic Press, London and NewYork. 1982: 123–139

LSQR  AN ALGORITHM FOR SPARSE LINEAREQUATIONS AND SPARSE LEASTSQUARES
ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE
1982; 8 (1): 4371
View details for Web of Science ID A1982NH42200005

A NOTE ON A SUFFICIENTDECREASE CRITERION FOR A NONDERIVATIVE STEPLENGTH PROCEDURE
MATHEMATICAL PROGRAMMING
1982; 23 (3): 349352
View details for Web of Science ID A1982NX42400006

A PROJECTED LAGRANGIAN ALGORITHM AND ITS IMPLEMENTATION FOR SPARSE NONLINEAR CONSTRAINTS
MATHEMATICAL PROGRAMMING STUDY
1982; 16 (MAR): 84117
View details for Web of Science ID A1982NR90100006

TOWARDS A GENERALIZED SINGULAR VALUE DECOMPOSITION
SIAM JOURNAL ON NUMERICAL ANALYSIS
1981; 18 (3): 398405
View details for Web of Science ID A1981LT69800003
 QPbased methods for largescale nonlinearly constrained optimization Nonlinear Programming 4 edited by Mangasarian, O., L., Meyer, R., R., M., S. Academic Press London and New York. 1981: 57–98
 A numerical investigation of ellipsoid algorithms for largescale linear programming Largescale Linear Programming edited by Dantzig, G., B., Dempster, M., A.H., Kallio, M. axenburg, Austria. 1981: 487–509

ASPECTS OF MATHEMATICALMODELING RELATED TO OPTIMIZATION
APPLIED MATHEMATICAL MODELLING
1981; 5 (2): 7183
View details for Web of Science ID A1981LJ58400002
 Methods for largescale nonlinear optimization Electric PowerProblems: The Mathematical Challenge edited by Erisman, A., M., Neves, K., W., Dwarakanath, M., H. SIAM, Philadelphia. 1980: 352–377
 Sparse least squares by conjugate gradients: a comparison of preconditioning methods 1979

LARGESCALE LINEARLY CONSTRAINED OPTIMIZATION
MATHEMATICAL PROGRAMMING
1978; 14 (1): 4172
View details for Web of Science ID A1978EM23400004

LEASTSQUARES ESTIMATION OF DISCRETE LINEAR DYNAMICSYSTEMS USING ORTHOGONAL TRANSFORMATIONS
SIAM JOURNAL ON NUMERICAL ANALYSIS
1977; 14 (2): 180193
View details for Web of Science ID A1977DC77900002

NONLINEAR OPTIMIZATION SUBJECT TO LINEARPROGRAMMING CONSTRAINTS
SIAM PUBLICATIONS. 1976: 825–26
View details for Web of Science ID A1976CH15400157
 A fast, stable implementation of the simplex method using BartelsGolub updating Sparse Matrix Computations edited by Bunch, J., R., Rose, D., J. Academic Press. 1976: 213–226
 The complexity of LU updating in the simplex method The Complexity of Computational Problem Solving edited by Brent, R., P. University of Queensland Press. 1976: 214–230

SOLUTION OF SPARSE INDEFINITE SYSTEMS OF LINEAR EQUATIONS
SIAM JOURNAL ON NUMERICAL ANALYSIS
1975; 12 (4): 617629
View details for Web of Science ID A1975AN11000008

METHODS FOR COMPUTING AND MODIFYING LDV FACTORS OF A MATRIX
MATHEMATICS OF COMPUTATION
1975; 29 (132): 10511077
View details for Web of Science ID A1975AU76800010
 Methods for computing and modifying the LDV factors of a matrix Math. Comput. 1975; 29: 10511077

METHODS FOR MODIFYING MATRIX FACTORIZATIONS
MATHEMATICS OF COMPUTATION
1974; 28 (126): 505535
View details for Web of Science ID A1974T209600011
 Numerical stability in largescale linear programming Approximation and Accuracy edited by deHoog, F., R., Jarvis, C., L. University of Queensland Press. 1973: 144–158
 Descent methods for minimization Optimization edited by Ryan, D., M. University of Queensland Press. 1972: 221–237
 Linear least squares and quadratic programming Integer and Nonlinear Programming edited by Abadie, J. NorthHolland, Amsterdam. 1970: 229–256
 Numerical techniques in mathematical programming Nonlinear Programming edited by Rosen, J., B., Mangasarian, O., L., Ritter, K. Academic Press, London and New York. 1970: 123–176
 Sparse least squares by conjugate gradients: a comparison of preconditioning methods edited by J. University of Waterloo, Waterloo, Ontario, Canada