Yuji Nakatsukasa
Current Organizations / Institutes
The organizations/institutes from the scholar's latest article
Previous Organizations / Institutes
The organizations / institutes from the scholar's other articles in the past ten years
Research Subjects
Yuji Nakatsukasa has a total of 113 co-authors. Here are the top co-authors by number of shared publications:
- Lloyd N. Trefethen
School of Engineering and Applied Sciences, Harvard University, Cambridge, USA
8 shared publications - Vanni Noferini
Department of Mathematics and Systems Analysis Aalto University P.O. Box 11100 FI-00076 Aalto Finland
8 shared publications - Alex Townsend
Department of Mathematics, Cornell University, Ithaca, USA
6 shared publications Lloyd N. Trefethen
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138 USA
6 shared publications- Paul Goulart
Department of Engineering Science University of Oxford OX1 3PJ Oxford UK
6 shared publications - Taejun Park
Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK
6 shared publications
Latest Publications
- Open Access
- 19 June 2026
- Published by Cambridge University Press (CUP) in Acta Numerica
- Vol. 35, 459-601
Applications of AAA rational approximation
The AAA algorithm for rational approximation is employed to illustrate applications of rational functions all across numerical analysis. For example, rational functions enable strikingly effective methods for numerical differentiation and integration; interpolation of equispaced data; locating zeros, poles, and branch points; analytic continuation; computation of inverse functions; imputation of missing data; computation of nonlinear eigenvalues and resonances; model order reduction/reduced order modelling; solution of Wiener–Hopf and Hilbert transform problems; conformal mapping; and solution of the two-dimensional Laplace, biharmonic and Helmholtz equations. The paper also surveys various improvements and generalizations of the AAA algorithm that have developed since its original appearance in 2018.
- 19 June 2026
- Published by Cambridge University Press (CUP) in Acta Numerica
- Vol. 35, 459-601
Open AccessApplications of AAA rational approximation
The AAA algorithm for rational approximation is employed to illustrate applications of rational functions all across numerical analysis. For example, rational functions enable strikingly effective methods for numerical differentiation and integration; interpolation of equispaced data; locating zeros, poles, and branch points; analytic continuation; computation of inverse functions; imputation of missing data; computation of nonlinear eigenvalues and resonances; model order reduction/reduced order modelling; solution of Wiener–Hopf and Hilbert transform problems; conformal mapping; and solution of the two-dimensional Laplace, biharmonic and Helmholtz equations. The paper also surveys various improvements and generalizations of the AAA algorithm that have developed since its original appearance in 2018.
- Open Access
- 27 February 2026
- Published by Elsevier in Applied and Computational Harmonic Analysis
- Vol. 83
Approximating sparse matrices and their functions using matrix-vector products
- 27 February 2026
- Published by Elsevier in Applied and Computational Harmonic Analysis
- Vol. 83
Open AccessApproximating sparse matrices and their functions using matrix-vector products
- 23 February 2026
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Matrix Analysis and Applications
- Vol. 47 (1) , 308-329
Sharp Error Bounds for Approximate Eigenvalues and Singular Values from Subspace Methods
- 23 February 2026
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Matrix Analysis and Applications
- Vol. 47 (1) , 308-329
Sharp Error Bounds for Approximate Eigenvalues and Singular Values from Subspace Methods
- Open Access
Fast, High-Accuracy, Randomized Nullspace Computations for Tall Matrices
In this paper, we develop RLOBPCG, an efficient method for computing a small number of singular triplets corresponding to the smallest singular values of large, tall matrices. The algorithm combines randomized preconditioner from the sketch-and-precondition techniques with the LOBPCG eigensolver: a small sketch is used to construct a high-quality preconditioner, and LOBPCG is run on the Gram matrix to refine the singular vector. Under the standard subspace embedding assumption and a modest singular value gap between the two smallest singular values, we prove that RLOBPCG converges geometrically to the minimum singular vector. In numerical experiments, RLOBPCG achieves near-optimal accuracy on matrices with up to $10^6$ rows, outperforming classical LOBPCG and Lanczos methods by a speedup of up to $12\times$ and maintaining robustness when other iterative methods fail to converge.Open AccessFast, High-Accuracy, Randomized Nullspace Computations for Tall Matrices
In this paper, we develop RLOBPCG, an efficient method for computing a small number of singular triplets corresponding to the smallest singular values of large, tall matrices. The algorithm combines randomized preconditioner from the sketch-and-precondition techniques with the LOBPCG eigensolver: a small sketch is used to construct a high-quality preconditioner, and LOBPCG is run on the Gram matrix to refine the singular vector. Under the standard subspace embedding assumption and a modest singular value gap between the two smallest singular values, we prove that RLOBPCG converges geometrically to the minimum singular vector. In numerical experiments, RLOBPCG achieves near-optimal accuracy on matrices with up to $10^6$ rows, outperforming classical LOBPCG and Lanczos methods by a speedup of up to $12\times$ and maintaining robustness when other iterative methods fail to converge. - 5 December 2025
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Matrix Analysis and Applications
- Vol. 46 (4) , 2614-2634
Matrix Perturbation Analysis of Methods for Extracting Singular Values from Approximate Singular Subspaces
Given (orthonormal) approximations [Formula: see text] and [Formula: see text] to the left and right subspaces spanned by the leading singular vectors of a matrix [Formula: see text], we discuss methods to approximate the leading singular values of [Formula: see text] and study their accuracy. In particular, we focus our analysis on the generalized Nyström approximation, as surprisingly, it is able to obtain significantly better accuracy than classical methods, namely Rayleigh–Ritz and (one-sided) projected SVD. A key idea of the analysis is to view the methods as finding the exact singular values of a perturbation of [Formula: see text]. In this context, we derive a matrix perturbation result that exploits the structure of such [Formula: see text] block matrix perturbation. This leads to sharp bounds that predict well the approximation error trends and explain the difference in the behavior of these methods. Finally, we present an approach to derive an a posteriori version of those bounds, which are more amenable to computation in practice.
- 5 December 2025
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Matrix Analysis and Applications
- Vol. 46 (4) , 2614-2634
Matrix Perturbation Analysis of Methods for Extracting Singular Values from Approximate Singular Subspaces
Given (orthonormal) approximations [Formula: see text] and [Formula: see text] to the left and right subspaces spanned by the leading singular vectors of a matrix [Formula: see text], we discuss methods to approximate the leading singular values of [Formula: see text] and study their accuracy. In particular, we focus our analysis on the generalized Nyström approximation, as surprisingly, it is able to obtain significantly better accuracy than classical methods, namely Rayleigh–Ritz and (one-sided) projected SVD. A key idea of the analysis is to view the methods as finding the exact singular values of a perturbation of [Formula: see text]. In this context, we derive a matrix perturbation result that exploits the structure of such [Formula: see text] block matrix perturbation. This leads to sharp bounds that predict well the approximation error trends and explain the difference in the behavior of these methods. Finally, we present an approach to derive an a posteriori version of those bounds, which are more amenable to computation in practice.
- Open Access
SubApSnap: Solving parameter-dependent linear systems with a snapshot and subsampling
A growing number of problems in computational mathematics can be reduced to the solution of many linear systems that are related, often depending smoothly or slowly on a parameter $p$, that is, $A(p)x(p)=b(p)$. We introduce an efficient algorithm for solving such parameter-dependent linear systems for many values of $p$. The algorithm, which we call SubApSnap (for \emph{Sub}sampled $A(p)$ times \emph{Snap}shot), is based on combining ideas from model order reduction and randomised linear algebra: namely, taking a snapshot matrix, and solving the resulting tall-skinny least-squares problems using a subsampling-based dimension-reduction approach. We show that SubApSnap is a strict generalisation of the popular DEIM algorithm in nonlinear model order reduction. SubApSnap is a sublinear-time algorithm, and once the snapshot and subsampling are determined, it solves $A(p_*)x(p_*)=b(p_*)$ for a new value of $p_*$ at a dramatically improved speed: it does not even need to read the whole matrix $A(p_*)$ to solve the linear system for a new value of $p_*$. We prove under natural assumptions that, given a good subsampling and snapshot, SubApSnap yields solutions with small residual for all parameter values of interest. We illustrate the efficiency and performance of the algorithm with problems arising in PDEs, model reduction, and kernel ridge regression, where SubApSnap achieves speedups of many orders of magnitude over a standard solution; for example over $20,000\times$ for a $10^7\times 10^7$ problem, while providing good accuracy.Open AccessSubApSnap: Solving parameter-dependent linear systems with a snapshot and subsampling
A growing number of problems in computational mathematics can be reduced to the solution of many linear systems that are related, often depending smoothly or slowly on a parameter $p$, that is, $A(p)x(p)=b(p)$. We introduce an efficient algorithm for solving such parameter-dependent linear systems for many values of $p$. The algorithm, which we call SubApSnap (for \emph{Sub}sampled $A(p)$ times \emph{Snap}shot), is based on combining ideas from model order reduction and randomised linear algebra: namely, taking a snapshot matrix, and solving the resulting tall-skinny least-squares problems using a subsampling-based dimension-reduction approach. We show that SubApSnap is a strict generalisation of the popular DEIM algorithm in nonlinear model order reduction. SubApSnap is a sublinear-time algorithm, and once the snapshot and subsampling are determined, it solves $A(p_*)x(p_*)=b(p_*)$ for a new value of $p_*$ at a dramatically improved speed: it does not even need to read the whole matrix $A(p_*)$ to solve the linear system for a new value of $p_*$. We prove under natural assumptions that, given a good subsampling and snapshot, SubApSnap yields solutions with small residual for all parameter values of interest. We illustrate the efficiency and performance of the algorithm with problems arising in PDEs, model reduction, and kernel ridge regression, where SubApSnap achieves speedups of many orders of magnitude over a standard solution; for example over $20,000\times$ for a $10^7\times 10^7$ problem, while providing good accuracy. - Open Access
Instability of the Sherman-Morrison formula and stabilization by iterative refinement
Owing to its simplicity and efficiency, the Sherman-Morrison (SM) formula has seen widespread use across various scientific and engineering applications for solving rank-one perturbed linear systems of the form $(A+uv^T)x = b$. Although the formula dates back at least to 1944, its numerical stability properties have remained an open question and continue to be a topic of current research. We analyze the backward stability of the SM, demonstrate its instability in a scenario increasingly common in scientific computing and address an open question posed by Nick Higham on the proportionality of the backward error bound to the condition number of $A$. We then incorporate fixed-precision iterative refinement into the SM framework reusing the previously computed decompositions and prove that, under reasonable assumptions, it achieves backward stability without sacrificing the efficiency of the SM formula. While our theory does not prove the SM formula with iterative refinement always outputs a backward stable solution, empirically it is observed to eventually produce a backward stable solution in all our numerical experiments. We conjecture that with iterative refinement, the SM formula yields a backward stable solution provided that $κ_2(A), κ_2(A+uv^T)$ are both bounded safely away from $ε_M^{-1}$, where $ε_M$ is the unit roundoff.Open AccessInstability of the Sherman-Morrison formula and stabilization by iterative refinement
Owing to its simplicity and efficiency, the Sherman-Morrison (SM) formula has seen widespread use across various scientific and engineering applications for solving rank-one perturbed linear systems of the form $(A+uv^T)x = b$. Although the formula dates back at least to 1944, its numerical stability properties have remained an open question and continue to be a topic of current research. We analyze the backward stability of the SM, demonstrate its instability in a scenario increasingly common in scientific computing and address an open question posed by Nick Higham on the proportionality of the backward error bound to the condition number of $A$. We then incorporate fixed-precision iterative refinement into the SM framework reusing the previously computed decompositions and prove that, under reasonable assumptions, it achieves backward stability without sacrificing the efficiency of the SM formula. While our theory does not prove the SM formula with iterative refinement always outputs a backward stable solution, empirically it is observed to eventually produce a backward stable solution in all our numerical experiments. We conjecture that with iterative refinement, the SM formula yields a backward stable solution provided that $κ_2(A), κ_2(A+uv^T)$ are both bounded safely away from $ε_M^{-1}$, where $ε_M$ is the unit roundoff. - Open Access
Efficient Function Approximation Under Heteroskedastic Noise
Approximating a function $f(x)$ on $[-1,1]$ based on $N+1$ samples is a classical problem in numerical analysis. If the samples come with heteroskedastic noise depending on $x$ of variance $σ(x)^2$, an $O(N\log N)$ algorithm for this problem has not yet been found in the current literature. In this paper, we propose a method called HeteroChebtrunc, adapted from an algorithm named NoisyChebtrunc. Using techniques in high-dimensional probability, we show that with high probability, HeteroChebtrunc achieves a tighter infinity-norm error bound than NoisyChebtrunc under heteroskedastic noise. This algorithm runs in $O(N+\hat{N}\log \hat{N})$ operations, where $\hat{N}\ll N$ is a chosen parameter. While investigating the properties of HeteroChebtrunc, we also derive a high-probability non-asymptotic relative error bound on the sample variance estimator for subgaussian variables, which is potentially another result of broader interest. We provide numerical experiments to demonstrate the improved uniform error of our algorithm.Open AccessEfficient Function Approximation Under Heteroskedastic Noise
Approximating a function $f(x)$ on $[-1,1]$ based on $N+1$ samples is a classical problem in numerical analysis. If the samples come with heteroskedastic noise depending on $x$ of variance $σ(x)^2$, an $O(N\log N)$ algorithm for this problem has not yet been found in the current literature. In this paper, we propose a method called HeteroChebtrunc, adapted from an algorithm named NoisyChebtrunc. Using techniques in high-dimensional probability, we show that with high probability, HeteroChebtrunc achieves a tighter infinity-norm error bound than NoisyChebtrunc under heteroskedastic noise. This algorithm runs in $O(N+\hat{N}\log \hat{N})$ operations, where $\hat{N}\ll N$ is a chosen parameter. While investigating the properties of HeteroChebtrunc, we also derive a high-probability non-asymptotic relative error bound on the sample variance estimator for subgaussian variables, which is potentially another result of broader interest. We provide numerical experiments to demonstrate the improved uniform error of our algorithm. - Open Access
- 9 July 2025
- Published by Springer Nature in Numerische Mathematik
- p. 1-27
Polynomial approximation of noisy functions
Approximating a univariate function on the interval $$[-1,1]$$ with a polynomial is among the most classical problems in numerical analysis. When the function evaluations come with noise, a least-squares fit is known to reduce the effect of noise as more samples are taken. The generic algorithm for the least-squares problem requires $$O(Nn^2)$$ operations, where $$N+1$$ is the number of sample points and n is the degree of the polynomial approximant. This algorithm is unstable when n is large, for example $$n\gg \sqrt{N}$$ for equispaced sample points. In this study, we blend numerical analysis and statistics to introduce a stable and fast $$O(N\log N)$$ algorithm called NoisyChebtrunc based on Chebyshev interpolation. It has the same error reduction effect as least-squares and the convergence is spectral until the error reaches $$O(\sigma \sqrt{{n}/{N}})$$ , where $$\sigma$$ is the noise level, after which the error continues to decrease at the Monte-Carlo $$O(1/\sqrt{N})$$ rate. To determine the polynomial degree, NoisyChebtrunc employs a statistical criterion, namely Mallows’ $$C_p$$ . We analyze NoisyChebtrunc in terms of the variance and concentration in the infinity norm to the underlying noiseless function. These results show that with high probability the infinity-norm error is bounded by a small constant times $$\sigma \sqrt{{n}/{N}}$$ , when the noise is independent and follows a subgaussian or subexponential distribution. We illustrate the performance of NoisyChebtrunc with numerical experiments.- 9 July 2025
- Published by Springer Nature in Numerische Mathematik
- p. 1-27
Open AccessPolynomial approximation of noisy functions
Approximating a univariate function on the interval $$[-1,1]$$ with a polynomial is among the most classical problems in numerical analysis. When the function evaluations come with noise, a least-squares fit is known to reduce the effect of noise as more samples are taken. The generic algorithm for the least-squares problem requires $$O(Nn^2)$$ operations, where $$N+1$$ is the number of sample points and n is the degree of the polynomial approximant. This algorithm is unstable when n is large, for example $$n\gg \sqrt{N}$$ for equispaced sample points. In this study, we blend numerical analysis and statistics to introduce a stable and fast $$O(N\log N)$$ algorithm called NoisyChebtrunc based on Chebyshev interpolation. It has the same error reduction effect as least-squares and the convergence is spectral until the error reaches $$O(\sigma \sqrt{{n}/{N}})$$ , where $$\sigma$$ is the noise level, after which the error continues to decrease at the Monte-Carlo $$O(1/\sqrt{N})$$ rate. To determine the polynomial degree, NoisyChebtrunc employs a statistical criterion, namely Mallows’ $$C_p$$ . We analyze NoisyChebtrunc in terms of the variance and concentration in the infinity norm to the underlying noiseless function. These results show that with high probability the infinity-norm error is bounded by a small constant times $$\sigma \sqrt{{n}/{N}}$$ , when the noise is independent and follows a subgaussian or subexponential distribution. We illustrate the performance of NoisyChebtrunc with numerical experiments. - 21 May 2025
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Scientific Computing
- Vol. 47 (3) , A1858-A1887
Low-Rank Approximation of Parameter-Dependent Matrices via CUR Decomposition
A low-rank approximation of a parameter-dependent matrix [Formula: see text] is an important task in the computational sciences appearing, for example, in dynamical systems and compression of a series of images. In this work, we introduce AdaCUR, an efficient algorithm for computing a low-rank approximation of parameter-dependent matrices via CUR decomposition. The key idea for this algorithm is that for nearby parameter values, the column and row indices for the CUR decomposition can often be reused. AdaCUR is rank-adaptive, provides error control, and has complexity that compares favorably against existing methods. A faster algorithm, FastAdaCUR, is also given that prioritizes speed over accuracy, is rank-adaptive, and has complexity that is at most linear in the number of rows or columns, but without error control.
- 21 May 2025
- Published by Society for Industrial & Applied Mathematics (SIAM) in SIAM Journal on Scientific Computing
- Vol. 47 (3) , A1858-A1887
Low-Rank Approximation of Parameter-Dependent Matrices via CUR Decomposition
A low-rank approximation of a parameter-dependent matrix [Formula: see text] is an important task in the computational sciences appearing, for example, in dynamical systems and compression of a series of images. In this work, we introduce AdaCUR, an efficient algorithm for computing a low-rank approximation of parameter-dependent matrices via CUR decomposition. The key idea for this algorithm is that for nearby parameter values, the column and row indices for the CUR decomposition can often be reused. AdaCUR is rank-adaptive, provides error control, and has complexity that compares favorably against existing methods. A faster algorithm, FastAdaCUR, is also given that prioritizes speed over accuracy, is rank-adaptive, and has complexity that is at most linear in the number of rows or columns, but without error control.