Pith. sign in

REVIEW 3 major objections 4 minor 57 references

Factorization-free Orthogonal Projection onto the Positive Semidefinite Cone with Composite Polynomial Filtering

T0 review · 3 major / 4 minor · reviewed 2026-08-06 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Composite polynomial filtering can replace eigenvalue decomposition for PSD cone projection, hitting 1e-3 relative error in half precision with 22 matrix multiplications and about 10x speedup.

desk verdict A well-engineered GPU method for approximate PSD projection with a real 10x speedup at moderate accuracy, but the abstract overstates consistency and the spectral-norm rescaling bound rests on an unproven heuristic. read the letter →

arxiv 2507.09165 v2 pith:VHMS53JW submitted 2025-07-12 math.OC

Shucheng Kang , Haoyu Han , Antoine Groudiev , Heng Yang This is my paper · ORCID
classification math.OC MSC 65F3090C2265F1565Y05
keywords PSDconeprojectioncompositepolynomialfilteringmatrixsignfunctionminimaxapproximationlow-precisionarithmeticGPUaccelerationsemidefiniteprogrammingmatrix-matrixmultiplication
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper tries to establish that orthogonally projecting a large symmetric matrix onto the positive semidefinite (PSD) cone—the operation that zeroes negative eigenvalues—can be replaced by a cascade of low-degree polynomial filters evaluated only with matrix-matrix multiplications, with no eigenvalue decomposition at all. Such a GEMM-only implementation suits GPU tensor cores and low-precision arithmetic, and the paper reports a consistent relative error near $10^{-3}$ in half precision using 22 matrix multiplications, roughly $10\times$ faster than the standard cuSOLVER eigendecomposition on matrices of size 5000 to 20000. If true, this makes repeated PSD projections in first-order semidefinite programming substantially cheaper, since first-order solvers only need moderate accuracy (KKT residuals around $10^{-2}$ to $10^{-4}$). Embedding the filter as a warm-start in an ADMM solver, the authors find low-precision projections reach the same final accuracy while cutting projection time by up to an order of magnitude.

What carries the argument

The load-bearing object is the composite polynomial filter of Eq. (6), $p=f_T\circ f_{T-1}\circ\cdots\circ f_1$, with each odd component of degree at most 5 and with coefficients fixed offline. The paper evaluates $p$ by applying the components sequentially to the matrix, so a degree-5 cascade of depth 7 or 10 costs 22 or 31 GEMMs rather than $5^T$; this is what makes a high-degree filter cheap on tensor cores. Two design stages produce the coefficients: Algorithm 1, a sequential Remez procedure whose output is minimax optimal for the sign function (Theorem 1), and a gradient-descent refinement on the ReLU objective (15) initialized from those coefficients. The runtime pipeline also rescales the input by a spectral-norm upper bound from the Residual Norm Theorem (Theorem 2), obtained from a 20-step Lanczos run, and rescales intermediate iterates to tame half-precision instability.

What would settle it

Fix a symmetric matrix with two nearly equal dominant eigenvalues (say $A$ with eigenvalues $1$ and $0.9999$ plus a tail near $0$), run the algorithm's 20-step Lanczos estimate on $A^2$, and compute whether (18) holds with the resulting $(\sigma,q)$; if ever $\tilde\lambda < \|A\|_2$, reconstruct the projection and check whether the relative error jumps well above $10^{-3}$, which would confirm the rescaling heuristic fails on that input.

Watch

Extended reading notes

Core claim

The central discovery is that the spectral operator $\Pi_{\mathbb{S}_+^n}(X)=Q\max(\lambda,0)Q^T$, whose scalar generator is the ReLU function $f_{\mathrm{ReLU}}(x)=\max\{x,0\}$, can be approximated by a composite polynomial $\tilde f(x)=\frac12 x(1+\tilde f_T\circ\cdots\circ\tilde f_1(x))$, where each $\tilde f_t$ is a degree-5 odd polynomial. Because the coefficients are kept in composite form, evaluating the filter on a dense matrix costs only the sum of the component degrees in general matrix-matrix multiplications, not the product. The coefficients are obtained in two stages: a sequential Remez construction that is minimax optimal for the matrix sign function, followed by gradient-based refinement that directly minimizes the worst-case ReLU error. With a Lanczos-based spectral-norm bound to rescale the matrix before filtering, the resulting half-precision filter attains median relative errors of $4.86\times10^{-4}$, $3.70\times10^{-4}$, and $5.14\times10^{-4}$ for $n=5000$, $10000$, and $20000$, respectively, on B200 GPUs—about $10\times$ faster than cuSOLVER's double-precision eigenvalue decomposition.

Load-bearing premise

The paper's theoretical spectral-norm bound is only proven when the shift $\sigma$ is the eigenvalue of $A^2$ nearest $\lambda_{\max}$, yet the algorithm takes $\sigma$ from a 20-step Lanczos run and asserts without proof that this heuristic gives a valid upper bound in all experiments; if it ever fails, the spectrum is not rescaled into $[-1,1]$ and the polynomial approximation guarantee collapses.

Editorial extensions

If this is right

  • Large dense PSD cone projections become memory-bandwidth-bound problems solved by GEMMs, so they can tap tensor cores; at about $10^{-3}$ relative error in half precision, a $10000\times10000$ projection takes about 55 ms on a B200 GPU.
  • First-order SDP solvers can safely use low-precision projections, since a KKT residual target of $10^{-2}$ to $10^{-4}$ does not require 64-bit projection accuracy; warm-starting ADMM with the filter preserves final accuracy.
  • The same sign-function cascade transfers to any spectral operator built from $\mathrm{sign}$, including $|X|$ and matrix polar factors, so the construction is not limited to the PSD cone.
  • Because the method needs no factorization, it avoids the orthogonality requirements that make eigendecomposition hard to implement in low precision.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • A rigorous check of the Lanczos heuristic behind Theorem 2 (that the largest Ritz value of $A^2$ is closest to $\lambda_{\max}$) would turn an empirical rescaling rule into a deterministic guarantee; until then the bound in (18) is the component most likely to fail silently.
  • The paper's noted failure mode—dominant extremal eigenvalues with many small nonzero ones, as in the triw dataset—suggests a natural hybrid: use the polynomial filter on the bulk spectrum and a deflated low-rank correction for the extremal eigenpairs, which could cut the GEMM count below 22.
  • On hardware without tensor cores or with slower GEMM, the 10x advantage may shrink or reverse, since the comparison is specifically against GPU eigendecomposition; a Roofline-style analysis with the machine's GEMM-to-memory balance point would predict where the crossover lies.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper proposes a factorization-free method for projecting a symmetric matrix onto the positive semidefinite cone. The projection operator f_ReLU(x)=max{x,0} is approximated by a composite polynomial p=f_T∘...∘f_1, where each f_i is a low-degree odd polynomial evaluated entirely through matrix-matrix multiplications. Coefficients are obtained in two stages: a minimax construction for the sign function following Lee et al. [27], followed by gradient-based refinement directly targeting f_ReLU. The authors implement the filter in CUDA with FP16 and FP32/BF16x9 emulation, benchmark it on 33 Matrix Depot datasets at sizes 5000, 10000, and 20000 on B200 and H100 GPUs, and report roughly 10x speedups over cuSOLVER FP64 with median FP16 relative errors around 5e-4. They also integrate the filter as a warm-start projection inside an ADMM solver for semidefinite programming. The central theoretical tool is a Lanczos-based spectral norm rescaling (Theorem 2) that is used to place the spectrum of the input in [-1,1] before applying the polynomial filter.

Significance. The paper's strengths are substantial: the implementations of both the filter and the benchmarks are open-sourced, the experimental comparison is extensive (33 datasets, three matrix sizes, five runs, two GPU architectures), and the idea of replacing eigenvalue decomposition with a short GEMM cascade is a useful and timely direction for low-accuracy, large-scale PSD projections. The SDP warm-start experiments give practical evidence that lower-accuracy projections can be useful in first-order solvers. However, the advertised 'consistent 10^-3' accuracy is not supported by the paper's own mean-error statistics, and the rescaling step on which all accuracy claims rest is only heuristic: Theorem 2's hypothesis is not verified by the algorithm. If the authors can certify or properly qualify the rescaling step and reconcile the abstract with the benchmark tables, the contribution would be a solid algorithms/systems paper.

major comments (3)
  1. [3.4, Theorem 2 and Algorithm 2] The rescaling step is load-bearing but not certified. Theorem 2 requires that sigma be the eigenvalue of A^2 closest to lambda_1 = lambda_max(A^2), but Algorithm 2 chooses sigma as the largest Ritz value of a 20-step Lanczos run and never checks this hypothesis. The residual inequality (19) bounds the distance from sigma to the eigenvalue nearest sigma, not to lambda_1; hence lambda tilde in (20) can under-estimate ||A||_2 if the Ritz value happens to be closer to a lower eigenvalue. Because the degree-5 component polynomials have coefficients as large as 22.6 in absolute value and grow rapidly outside [-1,1] (e.g., f_tilde_1(1.05) is about 2.7 from Table 2, and later components amplify further), even a modest under-estimate of the spectral norm can produce large projection errors. The statement in Section 3.4, item 1, that the heuristic 'yields a valid upper bound in all our extensive experiments' is empirical only and does not address this failure mode. Please provide a rigorous certificate for the scaling (for example, a guaranteed upper bound used as a fallback, or a proof that the chosen Lanczos pair satisfies the closest-eigenvalue condition), or explicitly present the scaling as a heuristic and adjust the theoretical claims accordingly.
  2. [Abstract; Tables 3-5] The abstract and Section 1 claim that Composite FP16 'consistently' attains a relative error of 10^-3, but the B200 mean relative errors in Tables 3, 4, and 5 are 9.53e-4, 3.45e-3, and 8.59e-3 for n=5000, 10000, and 20000, respectively. The medians (4.86e-4, 3.70e-4, 5.14e-4) do support a typical 10^-3 accuracy, but the mean error degrades by an order of magnitude at large n, driven by the type of outliers discussed in Section 4.1. Please revise the abstract and contribution statements to report median or typical accuracy, quantify the tail (e.g., maximum errors or worst-case datasets), or remove the word 'consistent'.
  3. [3.1, proof of Theorem 1] The proof of Theorem 1 as written has a reversed inequality. For any feasible polynomial sequence w, the optimal value E_min satisfies E_min <= max_x |w_T∘...∘w_1(x)-f_sign(x)|, so the displayed inequality E_min >= max_x |w_T∘...∘w_1(x)-f_sign(x)| is not a valid consequence and cannot lead to the claimed contradiction. Since the theorem is already established by [27, Theorem 2] and independently by [1, Theorem 4.1], please either correct the proof or replace it with a citation to those results.
minor comments (4)
  1. [Section 4.1, Figure 2] The legend of Figure 2 omits cuSOLVER FP64 because it is used as the ground truth for the error metric; please state this explicitly in the caption to avoid confusion.
  2. [Tables 3-5 and Figures 4-9] The distinction between 'Composite FP32' and 'Composite FP32 (em.)' is only explained in the text; please clarify in the table captions and figure legends that the emulated variant uses BF16x9 on B200 while the non-emulated variant uses native FP32 on H100.
  3. [Section 3.3, Tables 1 and 2] The column headers f_star_single and f_star_half are misleading because the formulas shown are for the individual component polynomials f_star_t, not for the full composite polynomial; please rename the columns or add a clarifying note.
  4. [Algorithm 2, line 7] The returned expression lambda_tilde * (1/2) X0 (I_n + X_T) is algebraically equal to (1/2) X (I_n + X_T) because X0 = X/lambda_tilde; please state this cancellation explicitly, as it clarifies how the scaling and unscaling interact with the sign approximation.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: coefficients are fitted to scalar approximation problems and validated against external cuSOLVER benchmarks, with no parameter tuned on the benchmark data.

full rationale

The derivation chain is not circular. The composite polynomial coefficients are produced by a two-stage offline procedure: Stage I minimizes the scalar minimax error for the sign function on [−1,1] via a sequential Remez algorithm, and Stage II locally refines these coefficients against the scalar ReLU objective on a dense grid of floating-point representable points. Neither stage uses matrix projection data; the matrix benchmarks (relative error vs. cuSOLVER FP64) are entirely external validation, so the reported 10^-3 accuracy is an independent empirical result rather than a fitted parameter renamed as a prediction. The minimax optimality of Stage I is imported from external references [27] and [1], not from the authors' own prior work, and the paper even includes a proof sketch for completeness. The spectral rescaling in Algorithm 2 uses a 20-step Lanczos Ritz pair as a heuristic choice of σ; the fact that the hypothesis of Theorem 2 (σ closest to λ_max(A²)) is not verified is a genuine correctness risk, but it is not circularity—it is an unproven assumption about a standard external residual-norm bound (Parlett), explicitly acknowledged in Section 3.4. Self-citations appear only in related-work and SDP-context discussion (e.g., [21], [23]) and are not load-bearing for the core factorization-free projection derivation. The final reconstruction X0(I+X_T)/2 is the standard algebraic identity ReLU(x) = x(1+sign(x))/2, not a definitional smuggling of the result. Overall, the paper is self-contained against external benchmarks and the central claims rest on measured performance, so no circular step is present.

Assumptions & free parameters 6 free parameters · 5 assumptions · 0 invented entities

The central claim rests on standard matrix spectral theory, on imported minimax results from [27] and [1], and on one heuristic (Lanczos-based scaling bound) that is not rigorously justified. The polynomial coefficients themselves are fitted to the approximation problem, which is legitimate but means the reported scalar errors are training errors.

free parameters (6)
  • Composite polynomial coefficients (single precision, T=10) = Listed in Table 1; e_float = 8.7023e-6
    Optimized by Stage I Remez for sign and Stage II gradient descent for ReLU on [-1,1]. These are the core design parameters of the filter.
  • Composite polynomial coefficients (half precision, T=7) = Listed in Table 2; e_float = 4.9233e-5
    Same optimization procedure as above, for half-precision arithmetic.
  • Epsilon (sign approximation gap) = 1e-3
    Hyperparameter in (9) that defines the interval [-1,-epsilon] U [epsilon,1] for sign approximation; chosen following [1,27].
  • Number of stages T and per-stage degree d_t = T=10 (single), T=7 (half); d_t=5 for all t
    Fixed computational budget in GEMMs; values guided by [1,27].
  • Post-iteration rescaling factors = 1/1.01 (FP16), 1/1.001 (FP32 first 8 iterations)
    Hand-tuned constants to stabilize low-precision evaluation of the composite polynomial.
  • Lanczos steps for spectral norm estimate = 20
    Choice of the number of Lanczos iterations to estimate the largest Ritz pair; running time is declared negligible.
assumptions (5)
  • standard math Spectral decomposition of symmetric matrices and the closed-form PSD projection (3) from Higham [17].
    Used as the definition of the target operator; standard linear algebra.
  • standard math For any polynomial p, p(X) = Q diag(p(lambda_i)) Q^T for symmetric X.
    Fundamental property of matrix polynomials under orthogonal diagonalization; used in Eq. (5).
  • standard math Theorem 1: the sequential Remez algorithm produces the minimax optimal composite for the sign function (9), as proven in [27, Theorem 2].
    Relied upon for Stage I initialization. The proof imports [27, Theorem 2] and uniqueness of Remez steps.
  • ad hoc to paper The 20-step Lanczos largest Ritz pair provides a valid upper bound in Theorem 2 for all inputs encountered.
    The authors state 'we have found this heuristic to yield a valid upper bound in all our extensive experiments' (Section 3.4.1), but provide no proof that the closest-eigenvalue condition holds. This is load-bearing for correct rescaling.
  • domain assumption The discretized error e_float over single-precision values correctly reflects the true sup-norm error.
    The paper uses e_float(f) as a proxy for e(f) in Stage II and in reported error tables; it is an approximation, not an exact computation.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Factorization-free Orthogonal Projection onto the Positive Semidefinite Cone with Composite Polynomial Filtering." pith.science (2026). https://pith.science/paper/VHMS53JW

@misc{pith2026250709165,
  author       = {Pith},
  title        = {Pith review of: Factorization-free Orthogonal Projection onto the Positive Semidefinite Cone with Composite Polynomial Filtering},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/VHMS53JW}},
  note         = {Machine review of arXiv:2507.09165}
}
abstract

We propose a factorization-free method for orthogonal projection onto the positive semidefinite (PSD) cone, leveraging composite polynomial filtering. Inspired by recent advances in homomorphic encryption, our approach approximates the PSD cone projection operator using a carefully optimized composite polynomial evaluated exclusively via matrix-matrix multiplications. This approach enables efficient GPU implementations with low-precision arithmetic, significantly outperforming the classical PSD cone projection using state-of-the-art GPU-based eigenvalue decomposition solvers. Specifically, our method achieves a consistent relative error of $10^{-3}$ in half-precision arithmetic with only 22 matrix-matrix multiplications, providing roughly a $10\times$ speed-up over NVIDIA's cuSOLVER routines on various large-scale matrices. In single-precision arithmetic with emulation on B200 GPUs, our approach maintains competitive accuracy while achieving up to a $2\times$ speed-up. Consequently, for a $10,000 \times 10,000$ dense symmetric matrix, our method requires approximately $55$ ms in half-precision and $400$ ms in single-precision arithmetic on B200 GPUs. Integration into a first-order semidefinite programming solver confirms that our low-precision projections reliably yield solutions of moderate accuracy.

Figures

Figures reproduced from arXiv: 2507.09165 by the authors.

Figure 1
Figure 1. Boxplots for different PSD cone projection methods’ [PITH_FULL_IMAGE:figures/full_fig_p011_1.png] view at source ↗
Figure 2
Figure 2. Boxplots for different PSD cone projection methods’ [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗
Figure 3
Figure 3. Numerical results for warm-starting ADMM with low-precision PSD-cone projection via composite [PITH_FULL_IMAGE:figures/full_fig_p014_3.png] view at source ↗
Figures from the paper (6 more)
Figure 4
Figure 4. Figure 4: Boxplots for different PSD cone projection methods’ [PITH_FULL_IMAGE:figures/full_fig_p016_4.png]
Figure 5
Figure 5. Figure 5: Boxplots for different PSD cone projection methods’ [PITH_FULL_IMAGE:figures/full_fig_p016_5.png]
Figure 6
Figure 6. Figure 6: PSD projection time on B200 GPU [PITH_FULL_IMAGE:figures/full_fig_p018_6.png]
Figure 7
Figure 7. Figure 7: PSD projection time on H100 GPU [PITH_FULL_IMAGE:figures/full_fig_p019_7.png]
Figure 8
Figure 8. Figure 8: PSD projection error on B200 GPU [PITH_FULL_IMAGE:figures/full_fig_p020_8.png]
Figure 9
Figure 9. Figure 9: PSD projection error on H100 GPU [PITH_FULL_IMAGE:figures/full_fig_p021_9.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

57 extracted references · 50 canonical work pages

  1. [27]

    Minimax approximation of sign function by composite polynomial for homomorphic comparison.IEEE Transactions on Dependable and Secure Computing, 19(6):3711–3727, 2022

    Eunsang Lee, Joon Woo Lee, Jong Seon No, and Young Sik Kim. Minimax approximation of sign function by composite polynomial for homomorphic comparison.IEEE Transactions on Dependable and Secure Computing, 19(6):3711–3727, 2022

    work page 2022

  2. [1]

    The polar express: Optimal matrix sign methods and their application to the muon algorithm.arXiv preprint arXiv:2505.16932, 2025

    Noah Amsel, David Persson, Christopher Musco, and Robert Gower. The polar express: Optimal matrix sign methods and their application to the muon algorithm.arXiv preprint arXiv:2505.16932, 2025

    arXiv 2025

  3. [2]

    SIAM, 1999

    Edward Anderson, Zhaojun Bai, Christian Bischof, L Susan Blackford, James Demmel, Jack Dongarra, Jeremy Du Croz, Anne Greenbaum, Sven Hammarling, Alan McKenney, et al.LAPACK users’ guide. SIAM, 1999

    work page 1999

  4. [3]

    Boosting projective methods for quantum process and detector tomography.Physical Review Research, 7:013208, 2025

    Júlia Barberà-Rodríguez, Leonardo Zambrano, Antonio Acín, and Donato Farina. Boosting projective methods for quantum process and detector tomography.Physical Review Research, 7:013208, 2025

    work page 2025

  5. [4]

    Squeezing 1-2% efficiency gains out of Muon by optimizing the Newton-Schulz coefficients.https://leloykun.github.io/ponder/ muon-opt-coeffs/, 2024

    Franz Louis Cesista, Jiacheng You, and Keller Jordan. Squeezing 1-2% efficiency gains out of Muon by optimizing the Newton-Schulz coefficients.https://leloykun.github.io/ponder/ muon-opt-coeffs/, 2024. Accessed: 2025-06-29

    work page 2024

  6. [5]

    A stable scaling of Newton-Schulz for improving the sign function compu- tation of a hermitian matrix

    Jie Chen and Edmond Chow. A stable scaling of Newton-Schulz for improving the sign function compu- tation of a hermitian matrix. Technical Report ANL/MCS-P5059-0114, Argonne National Laboratory, 2014

    work page 2014

  7. [6]

    Denman and Alex N

    Eugene D. Denman and Alex N. Beavers. The matrix sign function and computations in systems. Applied Mathematics and Computation, 2(1):63–94, 1976

    work page 1976

  8. [7]

    Spectral operators of matrices.Mathematical Programming, 168:509–531, 2018

    Chao Ding, Defeng Sun, Jie Sun, and Kim-Chuan Toh. Spectral operators of matrices.Mathematical Programming, 168:509–531, 2018

    work page 2018

Show all 57 references
  1. [8]

    Pelejo, Yuen-Lam Voronin, and Henry Wolkowicz

    Dmitriy Drusvyatskiy, Chi-Kwong Li, Diane C. Pelejo, Yuen-Lam Voronin, and Henry Wolkowicz. Projection methods for quantum channel construction.Quantum Information Processing, 14:3075–3096, 2015

    2015

  2. [9]

    Francisco and Douglas S

    Juliano B. Francisco and Douglas S. Gonçalves. A fixed-point method for approximate projection onto the positive semidefinite cone.Linear Algebra and its Applications, 523:59–78, 2017

    2017

  3. [10]

    Cosmo: A conic operator splitting method for convex conic problems.Journal of Optimization Theory and Applications, 190(3):779–810, 2021

    Michael Garstka, Mark Cannon, and Paul Goulart. Cosmo: A conic operator splitting method for convex conic problems.Journal of Optimization Theory and Applications, 190(3):779–810, 2021

    2021

  4. [11]

    Networks and the best approximation property.Biological Cyber- netics, 63(3):169–176, 1990

    Federico Girosi and Tomaso Poggio. Networks and the best approximation property.Biological Cyber- netics, 63(3):169–176, 1990

    1990

  5. [12]

    MAGMA: A breakthrough in solvers for eigenvalue problems

    Abdulla Haidar, Stanimire Tomov, Ichitaro Yamazaki, Romain Solca, Thomas Schulthess, Tingxing Dong, and Jack Dongarra. MAGMA: A breakthrough in solvers for eigenvalue problems. White paper, GPU Technology Conference (GTC’12), May 2012

    2012

  6. [13]

    Building rome with convex optimization

    Haoyu Han and Heng Yang. Building rome with convex optimization. InRobotics: Science and Systems (RSS), 2025

    2025

  7. [14]

    Projection methods for conic feasibility problems: applications to polynomial sum-of-squares decompositions.Optimization Methods & Software, 26(1):23–46, 2011

    Didier Henrion and Jérôme Malick. Projection methods for conic feasibility problems: applications to polynomial sum-of-squares decompositions.Optimization Methods & Software, 26(1):23–46, 2011

    2011

  8. [15]

    Projection methods in conic optimization

    Didier Henrion and Jérôme Malick. Projection methods in conic optimization. In Miguel F. Anjos and Jean B. Lasserre, editors,Handbook on Semidefinite, Conic and Polynomial Optimization, volume 166 ofInternational Series in Operations Research & Management Science. Springer, 2012

    2012

  9. [16]

    Nicholas J. Higham. Computing the polar decomposition—with applications.SIAM Journal on Scien- tific and Statistical Computing, 7(4):1160–1174, 1986

    1986

  10. [17]

    Computing a nearest symmetric positive semidefinite matrix.Linear algebra and its applications, 103:103–118, 1988

    Nicholas J Higham. Computing a nearest symmetric positive semidefinite matrix.Linear algebra and its applications, 103:103–118, 1988. 22

    1988

  11. [18]

    Nicholas J. Higham. Computing the nearest correlation matrix—a problem from finance.IMA Journal of Numerical Analysis, 22(3):329–343, 2002

    2002

  12. [19]

    SocietyforIndustrialandApplied Mathematics, Philadelphia, PA, 2008

    NicholasJ.Higham.Functions of Matrices: Theory and Computation. SocietyforIndustrialandApplied Mathematics, Philadelphia, PA, 2008

    2008

  13. [20]

    Approximate projections onto the positive semidefinite cone using randomization.arXiv preprint arXiv:2410.19208, 2024

    Morgan Jones and James Anderson. Approximate projections onto the positive semidefinite cone using randomization.arXiv preprint arXiv:2410.19208, 2024

    2024 arXiv

  14. [21]

    Local linear convergence of the alternating direction method of multipliers for semidefinite programming under strict complementarity.arXiv preprint arXiv:2503.20142, 2025

    Shucheng Kang, Xin Jiang, and Heng Yang. Local linear convergence of the alternating direction method of multipliers for semidefinite programming under strict complementarity.arXiv preprint arXiv:2503.20142, 2025

    2025 arXiv

  15. [22]

    Global contact-rich planning with sparsity-rich semidefi- nite relaxations

    Shucheng Kang, Guorui Liu, and Heng Yang. Global contact-rich planning with sparsity-rich semidefi- nite relaxations. InRobotics: Science and Systems (RSS), 2025

    2025

  16. [23]

    Fast and certifiable trajectory optimization

    Shucheng Kang, Xiaoyang Xu, Jay Sarva, Ling Liang, and Heng Yang. Fast and certifiable trajectory optimization. InInternational Workshop on the Algorithmic Foundations of Robotics (WAFR), 2024

    2024

  17. [24]

    State polynomials: positivity, optimization and nonlinear bell inequalities.Mathematical Programming, 207(1):645–691, 2024

    Igor Klep, Victor Magron, Jurij Volčič, and Jie Wang. State polynomials: positivity, optimization and nonlinear bell inequalities.Mathematical Programming, 207(1):645–691, 2024

    2024

  18. [25]

    Andrew V. Knyazev. Toward the optimal preconditioned eigensolver: Locally optimal block precondi- tioned conjugate gradient method.SIAM Journal on Scientific Computing, 23(2):517–541, 2001

    2001

  19. [26]

    Min-maxandrobustpolynomialoptimization.Journal of Global Optimization, 51(1):1– 10, 2011

    JeanBLasserre. Min-maxandrobustpolynomialoptimization.Journal of Global Optimization, 51(1):1– 10, 2011

    2011

  20. [28]

    World Scientific, 2023

    Victor Magron and Jie Wang.Sparse polynomial optimization: theory and practice. World Scientific, 2023

    2023

  21. [29]

    A dual approach to semidefinite least-squares problems.SIAM Journal on Matrix Analysis and Applications, 26(1):272–284, 2004

    Jérôme Malick. A dual approach to semidefinite least-squares problems.SIAM Journal on Matrix Analysis and Applications, 26(1):272–284, 2004

    2004

  22. [30]

    Regularizationmethodsforsemidefinite programming.SIAM Journal on Optimization, 20(1):336–356, 2009

    JérômeMalick, JanezPovh, FranzRendl, andAngelikaWiegele. Regularizationmethodsforsemidefinite programming.SIAM Journal on Optimization, 20(1):336–356, 2009

    2009

  23. [31]

    Several sdp-codes on sparse and other sdp problems, 2006

    Hans D Mittelmann. Several sdp-codes on sparse and other sdp problems, 2006

    2006

  24. [32]

    Optimizing Halley’s iteration for computing the matrix polar decomposition.SIAM Journal on Matrix Analysis and Applications, 31(5):2700–2720, 2010

    Yuji Nakatsukasa, Zhaojun Bai, and François Gygi. Optimizing Halley’s iteration for computing the matrix polar decomposition.SIAM Journal on Matrix Analysis and Applications, 31(5):2700–2720, 2010

    2010

  25. [33]

    Donald J. Newman. Rational approximation to|x|.The Michigan Mathematical Journal, 11(1):11–14, 1964

    1964

  26. [34]

    Whitepaper

    NVIDIA Corporation.NVIDIA A100 Tensor Core GPU Architecture, 2020. Whitepaper

    2020

  27. [35]

    NVIDIA Corporation. Boosting matrix multiplication speed and flex- ibility with nvidia cublas 12.9.https://developer.nvidia.com/blog/ boosting-matrix-multiplication-speed-and-flexibility-with-nvidia-cublas-12-9/, May

  28. [36]

    Version 12.x

    NVIDIA Corporation.NVIDIA cuSOLVER Library, 2024. Version 12.x. 23

    2024

  29. [37]

    NVIDIA Developer Documen- tation

    NVIDIA Corporation.cuBLAS Library User Guide, Version 12.9, 2025. NVIDIA Developer Documen- tation

    2025

  30. [38]

    SCS: Splitting conic solver, version 3.2.4.https://github.com/cvxgrp/scs, November 2023

    Brendan O’Donoghue, Eric Chu, Neal Parikh, and Stephen Boyd. SCS: Splitting conic solver, version 3.2.4.https://github.com/cvxgrp/scs, November 2023

    2023

  31. [39]

    SIAM, Philadelphia, PA, 1998

    Beresford N Parlett.The symmetric eigenvalue problem. SIAM, Philadelphia, PA, 1998

    1998

  32. [40]

    E. Ya. Remez.General Computational Methods of Chebyshev Approximation: The Problems with Linear Realization. Kiev: Izdatel’stvo Akademii Nauk Ukrainskoi SSR, 1962. In Russian; English translation by J. H. McCabe, 1964, U.S. Department of Commerce, Office of Technical Services

    1962

  33. [41]

    Efficient semidefinite programming with approx- imate ADMM.Journal of Optimization Theory and Applications, pages 1–29, 2022

    Nikitas Rontsis, Paul Goulart, and Yuji Nakatsukasa. Efficient semidefinite programming with approx- imate ADMM.Journal of Optimization Theory and Applications, pages 1–29, 2022

    2022

  34. [42]

    Exploiting low-rank structure in semidefinite programming by approximate operator splitting.Optimization, 71(1):117–144, 2022

    Mario Souto, Joaquim D Garcia, and Álvaro Veiga. Exploiting low-rank structure in semidefinite programming by approximate operator splitting.Optimization, 71(1):117–144, 2022

    2022

  35. [43]

    Towards dense linear algebra for hybrid GPU accelerated manycore systems.Parallel Computing, 36(5-6):232–240, 2010

    Stanimire Tomov, Jack Dongarra, and Marc Baboulin. Towards dense linear algebra for hybrid GPU accelerated manycore systems.Parallel Computing, 36(5-6):232–240, 2010

    2010

  36. [44]

    SIAM, 2019

    Lloyd N Trefethen.Approximation theory and approximation practice, extended edition. SIAM, 2019

    2019

  37. [45]

    Vasily Volkov and James W. Demmel. Benchmarking GPUs to tune dense linear algebra. InProceedings of the 2008 ACM/IEEE Conference on Supercomputing (SC’08), pages 1–11. IEEE, 2008

    2008

  38. [46]

    Certifying ground-state properties of many-body systems.Physical Review X, 14(3):031006, 2024

    Jie Wang, Jacopo Surace, Irénée Frérot, Benoît Legat, Marc-Olivier Renou, Victor Magron, and Antonio Acín. Certifying ground-state properties of many-body systems.Physical Review X, 14(3):031006, 2024

    2024

  39. [47]

    Weinberger and Lawrence K

    Kilian Q. Weinberger and Lawrence K. Saul. Distance metric learning for large margin nearest neighbor classification.Journal of Machine Learning Research, 10:207–244, 2009

    2009

  40. [48]

    Alternating direction augmented Lagrangian methods for semidefinite programming.Mathematical Programming Computation, 2(3):203–230, 2010

    Zaiwen Wen, Donald Goldfarb, and Wotao Yin. Alternating direction augmented Lagrangian methods for semidefinite programming.Mathematical Programming Computation, 2(3):203–230, 2010

    2010

  41. [49]

    Roofline: An insightful visual performance model for multicore architectures.Communications of the ACM, 52(4):65–76, 2009

    Samuel Williams, Andrew Waterman, and David Patterson. Roofline: An insightful visual performance model for multicore architectures.Communications of the ACM, 52(4):65–76, 2009

    2009

  42. [50]

    Heng Yang, Ling Liang, Luca Carlone, and Kim-Chuan Toh. An inexact projected gradient method with rounding and lifting by nonlinear programming for solving rank-one semidefinite relaxation of polynomial optimization.Mathematical Programming, 201(1):409–472, 2023

    2023

  43. [51]

    Rounding error analysis of mixed precision block Householder QR algorithms.SIAM Journal on Scientific Computing, 43(3):A1723–A1753, 2021

    L Minah Yang, Alyson Fox, and Geoffrey Sanders. Rounding error analysis of mixed precision block Householder QR algorithms.SIAM Journal on Scientific Computing, 43(3):A1723–A1753, 2021

    2021

  44. [52]

    Liuqin Yang, Defeng Sun, and Kim-Chuan Toh. SDPNAL+: A majorized semismooth Newton-CG aug- mented Lagrangian method for semidefinite programming with nonnegative constraints.Mathematical Programming Computation, 7(3):331–366, 2015

    2015

  45. [53]

    Structured channel covariance estimation from limited samples for large antenna arrays.EURASIP Journal on Wireless Communications and Networking, 2023(24):1–34, 2023

    Tianyu Yang, Mahdi Barzegar Khalilsarai, Saeid Haghighatshoar, and Giuseppe Caire. Structured channel covariance estimation from limited samples for large antenna arrays.EURASIP Journal on Wireless Communications and Networking, 2023(24):1–34, 2023

    2023

  46. [54]

    Matrix depot: an extensible test matrix collection for julia

    Weijian Zhang and Nicholas J Higham. Matrix depot: an extensible test matrix collection for julia. PeerJ Computer Science, 2:e58, 2016

    2016

  47. [55]

    Fast ADMM for homogeneous self-dual embedding of sparse SDPs.IFAC-PapersOnLine, 50(1):8411–8416, 2017

    Yang Zheng, Giovanni Fantuzzi, Antonis Papachristodoulou, Paul Goulart, and Andrew Wynn. Fast ADMM for homogeneous self-dual embedding of sparse SDPs.IFAC-PapersOnLine, 50(1):8411–8416, 2017. 24

    2017

  48. [56]

    An iterative algorithm for computing the best estimate of an orthogonal matrix.SIAM Journal on Numerical Analysis, 8(2):358–364, 1971

    Åke Björck and Charles Bowie. An iterative algorithm for computing the best estimate of an orthogonal matrix.SIAM Journal on Numerical Analysis, 8(2):358–364, 1971. 25

    1971

  49. [2024]

    Accessed: 2025-06-14

    2025

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.