Many problems in linear algebra -- such as those arising from non-Hermitian physics and differential equations -- can be solved on a quantum computer by processing eigenvalues of the non-normal input matrices. However, the existing Quantum Singular Value Transformation (QSVT) framework is ill-suited to this task, as eigenvalues and singular values are different in general. We present a Quantum EigenValue Transformation (QEVT) framework for applying arbitrary polynomial transformations on eigenvalues of block-encoded non-normal operators, and a related Quantum EigenValue Estimation (QEVE) algorithm for operators with real spectra. QEVT has query complexity to the block encoding nearly recovering that of the QSVT for a Hermitian input, and QEVE achieves the Heisenberg-limited scaling for diagonalizable input matrices. As applications, we develop a linear differential equation solver with strictly linear time query complexity for average-case diagonalizable operators, as well as a ground state preparation algorithm that upgrades previous nearly optimal results for Hermitian Hamiltonians to diagonalizable matrices with real spectra. Underpinning our algorithms is an efficient method to prepare a quantum superposition of Faber polynomials, which generalize the nearly-best uniform approximation properties of Chebyshev polynomials to the complex plane. Of independent interest, we also develop techniques to generate $n$ Fourier coefficients with $\mathbf{O}(\mathrm{polylog}(n))$ gates compared to prior approaches with linear cost.
We introduce a family of identities that express general linear non-unitary evolution operators as a linear combination of unitary evolution operators, each solving a Hamiltonian simulation problem. This formulation can exponentially enhance the accuracy of the recently introduced linear combination of Hamiltonian simulation (LCHS) method [An, Liu, and Lin, Physical Review Letters, 2023]. For the first time, this approach enables quantum algorithms to solve linear differential equations with both optimal state preparation cost and near-optimal scaling in matrix queries on all parameters.
We present an efficient quantum circuit for block encoding pairing Hamiltonian often studied in nuclear physics. Our block encoding scheme does not require mapping the creation and annihilation operators to the Pauli operators and representing the Hamiltonian as a linear combination of unitaries. Instead, we show how to encode the Hamiltonian directly using controlled swap operations. We analyze the gate complexity of the block encoding circuit and show that it scales polynomially with respect to the number of qubits required to represent a quantum state associated with the pairing Hamiltonian. We also show how the block encoding circuit can be combined with the quantum singular value transformation to construct an efficient quantum circuit for approximating the density of states of a pairing Hamiltonian. The techniques presented can be extended to encode more general second-quantized Hamiltonians.
While non-Hermitian physics has attracted considerable attention, current studies are limited to small or classically solvable systems. Quantum computing, as a powerful eigensolver, have predominantly been applied to Hermitian domain, leaving their potential for studying non-Hermitian problems largely unexplored. We extend the power of quantum computing to general non-Hermitian eigenproblems. Our approach works for finding eigenvalues without extra constrains, or eigenvalues closest to specified points or lines, thus extending results for ground energy and energy gap problems for Hermitian matrices. Our algorithms have broad applications, and as examples, we consider two central problems in non-Hermitian physics. Firstly, our approach is the first to offer an efficient quantum solution to the witness of spontaneous $PT$-symmetry breaking, and provide provable, exponential quantum advantage. Secondly, our approach enables the estimation of Liouvillian gap, which is crucial for characterizing relaxation times. Our general approach can also find applications in many other areas, such as the study of Markovian stochastic processes. These results underscore the significance of our quantum algorithms for addressing practical eigenproblems across various disciplines.
Quantum computing has emerged as a promising avenue for achieving significant speedup, particularly in large-scale PDE simulations, compared to classical computing. One of the main quantum approaches involves utilizing Hamiltonian simulation, which is directly applicable only to Schrödinger-type equations. To address this limitation, Schrödingerisation techniques have been developed, employing the warped transformation to convert general linear PDEs into Schrödinger-type equations. However, despite the development of Schrödingerisation techniques, the explicit implementation of the corresponding quantum circuit for solving general PDEs remains to be designed. In this paper, we present detailed implementation of a quantum algorithm for general PDEs using Schrödingerisation techniques. We provide examples of the heat equation, and the advection equation approximated by the upwind scheme, to demonstrate the effectiveness of our approach. Complexity analysis is also carried out to demonstrate the quantum advantages of these algorithms in high dimensions over their classical counterparts.
High-dimensional fractional reaction-diffusion equations have numerous applications in the fields of biology, chemistry, and physics, and exhibit a range of rich phenomena. While classical algorithms have an exponential complexity in the spatial dimension, a quantum computer can produce a quantum state that encodes the solution with only polynomial complexity, provided that suitable input access is available. In this work, we investigate efficient quantum algorithms for linear and nonlinear fractional reaction-diffusion equations with periodic boundary conditions. For linear equations, we analyze and compare the complexity of various methods, including the second-order Trotter formula, time-marching method, and truncated Dyson series method. We also present a novel algorithm that combines the linear combination of Hamiltonian simulation technique with the interaction picture formalism, resulting in optimal scaling in the spatial dimension. For nonlinear equations, we employ the Carleman linearization method and propose a block-encoding version that is appropriate for the dense matrices that arise from the spatial discretization of fractional reaction-diffusion equations.
The well-conditioned multi-product formula (MPF), proposed by [Low, Kliuchnikov, and Wiebe, 2019], is a simple high-order time-independent Hamiltonian simulation algorithm that implements a linear combination of standard product formulas of low order. While the MPF aims to simultaneously exploit commutator scaling among Hamiltonians and achieve near-optimal time and precision dependence, its lack of a rigorous error bound on the nested commutators renders its practical advantage ambiguous. In this work, we conduct a rigorous complexity analysis of the well-conditioned MPF, demonstrating explicit commutator scaling and near-optimal time and precision dependence at the same time. Using our improved complexity analysis, we present several applications of practical interest where the MPF based on a second-order product formula can achieve a polynomial speedup in both system size and evolution time, as well as an exponential speedup in precision, compared to second-order and even higher-order product formulas. Compared to post-Trotter methods, the MPF based on a second-order product formula can achieve polynomially better scaling in system size, with only poly-logarithmic overhead in evolution time and precision.
Solving linear systems is at the foundation of many algorithms. Recently, quantum linear system algorithms (QLSAs) have attracted great attention since they converge to a solution exponentially faster than classical algorithms in terms of the problem dimension. However, low-complexity circuit implementations of the oracles assumed in these QLSAs constitute the major bottleneck for practical quantum speed-up in solving linear systems. In this work, we focus on the application of QLSAs for linear systems that are expressed as a low rank tensor sums, which arise in solving discretized PDEs. Previous works uses modified Krylov subspace methods to solve such linear systems with a per-iteration complexity being polylogarithmic of the dimension but with no guarantees on the total convergence cost. We propose a quantum algorithm based on the recent advances on adiabatic-inspired QLSA and perform a detailed analysis of the circuit depth of its implementation. We rigorously show that the total complexity of our implementation is polylogarithmic in the dimension, which is comparable to the per-iteration complexity of the classical heuristic methods.
In this article we introduce an algorithm for mitigating the adverse effects of noise on gradient descent in variational quantum algorithms. This is accomplished by computing a \emphregularized local classical approximation to the objective function at every gradient descent step. The computational overhead of our algorithm is entirely classical, i.e., the number of circuit evaluations is exactly the same as when carrying out gradient descent using the parameter-shift rules. We empirically demonstrate the advantages offered by our algorithm on randomized parametrized quantum circuits.
Hamiltonian simulation is a domain where quantum computers have the potential to outperform their classical counterparts due to their inherent quantum behavior. One of the main challenges of such quantum algorithms is up-scaling the system size, which is necessary to achieve meaningful quantum advantage. In this work, we present an approach to improve the scalability of eigenspace filtering for the ground state preparation of a given Hamiltonian. Our method aims to tackle limitations introduced by a small spectral gap and high degeneracy of low energy states. It is based on an adaptive sequence of eigenspace filtering through Quantum Eigenvalue Transformation of Unitary Matrices (QETU) followed by spectrum profiling. By combining our proposed algorithm with state-of-the-art phase estimation methods, we achieved good approximations for the ground state energy with local, two-qubit gate depolarizing probability up to $10^{-4}$. To demonstrate the key results in this work, we ran simulations with the transverse-field Ising Model on classical computers using Qiskit. We compare the performance of our approach with the static implementation of QETU and show that we can consistently achieve three to four orders of magnitude improvement in the absolute error rate.
The quantum Fourier transform (QFT), which can be viewed as a reindexing of the discrete Fourier transform (DFT), has been shown to be compressible as a low-rank matrix product operator (MPO) or quantized tensor train (QTT) operator. However, the original proof of this fact does not furnish a construction of the MPO with a guaranteed error bound. Meanwhile, the existing practical construction of this MPO, based on the compression of a quantum circuit, is not as efficient as possible. We present a simple closed-form construction of the QFT MPO using the interpolative decomposition, with guaranteed near-optimal compression error for a given rank. This construction can speed up the application of the QFT and the DFT, respectively, in quantum circuit simulations and QTT applications. We also connect our interpolative construction to the approximate quantum Fourier transform (AQFT) by demonstrating that the AQFT can be viewed as an MPO constructed using a different interpolation scheme.
Many promising quantum applications depend on the efficient quantum simulation of an exponentially large sparse Hamiltonian, a task known as sparse Hamiltonian simulation, which is fundamentally important in quantum computation. Although several theoretically appealing quantum algorithms have been proposed for this task, they typically require a black-box query model of the sparse Hamiltonian, rendering them impractical for near-term implementation on quantum devices. In this paper, we propose a technique named Hamiltonian embedding. This technique simulates a desired sparse Hamiltonian by embedding it into the evolution of a larger and more structured quantum system, allowing for more efficient simulation through hardware-efficient operations. We conduct a systematic study of this new technique and demonstrate significant savings in computational resources for implementing prominent quantum applications. As a result, we can now experimentally realize quantum walks on complicated graphs (e.g., binary trees, glued-tree graphs), quantum spatial search, and the simulation of real-space Schrödinger equations on current trapped-ion and neutral-atom platforms. Given the fundamental role of Hamiltonian evolution in the design of quantum algorithms, our technique markedly expands the horizon of implementable quantum advantages in the NISQ era.
This paper presents a quantum algorithm for the solution of prototypical second-order linear elliptic partial differential equations discretized by $d$-linear finite elements on Cartesian grids of a bounded $d$-dimensional domain. An essential step in the construction is a BPX preconditioner, which transforms the linear system into a sufficiently well-conditioned one, making it amenable to quantum computation. We provide a constructive proof demonstrating that, for any fixed dimension, our quantum algorithm can compute suitable functionals of the solution to a given tolerance $\mathtt{tol}$ with an optimal complexity of order $\mathtt{tol}^{-1}$ up to logarithmic terms, significantly improving over existing approaches. Notably, this approach does not rely on regularity of the solution and achieves quantum advantage over classical solvers in two dimensions, whereas prior quantum methods required at least four dimensions for asymptotic benefits. We further detail the design and implementation of a quantum circuit capable of executing our algorithm, present simulator results, and report numerical experiments on current quantum hardware, confirming the feasibility of preconditioned finite element methods for near-term quantum computing.
We analyze the Schrödingerisation method for quantum simulation of a general class of non-unitary dynamics with inhomogeneous source terms. The Schrödingerisation technique, introduced in \citeJLY22a,JLY23, transforms any linear ordinary and partial differential equations with non-unitary dynamics into a system under unitary dynamics via a warped phase transition that maps the equations into a higher dimension, making them suitable for quantum simulation. This technique can also be applied to these equations with inhomogeneous terms modeling source or forcing terms or boundary and interface conditions, and discrete dynamical systems such as iterative methods in numerical linear algebra, through extra equations in the system. Difficulty airses with the presense of inhomogeneous terms since it can change the stability of the original system. In this paper, we systematically study--both theoretically and numerically--the important issue of recovering the original variables from the Schrödingerized equations, even when the evolution operator contains unstable modes. We show that even with unstable modes, one can still construct a stable scheme, yet to recover the original variable one needs to use suitable data in the extended space. We analyze and compare both the discrete and continuous Fourier transforms used in the extended dimension, and derive corresponding error estimates, which allows one to use the more appropriate transform for specific equations. We also provide a smoother initialization for the Schrodödingerized system to gain higher order accuracy in the extended space. We homogenize the inhomogeneous terms with a stretch transformation, making it easier to recover the original variable. Our recovering technique also provides a simple and generic framework to solve general ill-posed problems in a computationally stable way.
This work proposes a protocol for Fermionic Hamiltonian learning. For the Hubbard model defined on a bounded-degree graph, the Heisenberg-limited scaling is achieved while allowing for state preparation and measurement errors. To achieve $\epsilon$-accurate estimation for all parameters, only $\tilde{\mathcal{O}}(\epsilon^{-1})$ total evolution time is needed, and the constant factor is independent of the system size. Moreover, our method only involves simple one or two-site Fermionic manipulations, which is desirable for experiment implementation.
Efficient error estimates for the Trotter product formula are central in quantum computing, mathematical physics, and numerical simulations. However, the Trotter error's dependency on the input state and its application to unbounded operators remains unclear. Here, we present a general theory for error estimation, including higher-order product formulas, with explicit input state dependency. Our approach overcomes two limitations of the existing operator-norm estimates in the literature. First, previous bounds are too pessimistic as they quantify the worst-case scenario. Second, previous bounds become trivial for unbounded operators and cannot be applied to a wide class of Trotter scenarios, including atomic and molecular Hamiltonians. Our method enables analytical treatment of Trotter errors in chemistry simulations, illustrated through a case study on the hydrogen atom. Our findings reveal: (i) for states with fat-tailed energy distribution, such as low-angular-momentum states of the hydrogen atom, the Trotter error scales worse than expected (sublinearly) in the number of Trotter steps; (ii) certain states do not admit an advantage in the scaling from higher-order Trotterization, and thus, the higher-order Trotter hierarchy breaks down for these states, including the hydrogen atom's ground state; (iii) the scaling of higher-order Trotter bounds might depend on the order of the Hamiltonians in the Trotter product for states with fat-tailed energy distribution. Physically, the enlarged Trotter error is caused by the atom's ionization due to the Trotter dynamics. Mathematically, we find that certain domain conditions are not satisfied by some states so higher moments of the potential and kinetic energies diverge. Our analytical error analysis agrees with numerical simulations, indicating that we can estimate the state-dependent Trotter error scaling genuinely.
Foundation models, such as large language models, have demonstrated success in addressing various language and image processing tasks. In this work, we introduce a multi-modal foundation model for scientific problems, named PROSE-PDE. Our model, designed for bi-modality to bi-modality learning, is a multi-operator learning approach which can predict future states of spatiotemporal systems while concurrently learning the underlying governing equations of the physical system. Specifically, we focus on multi-operator learning by training distinct one-dimensional time-dependent nonlinear constant coefficient partial differential equations, with potential applications to many physical applications including physics, geology, and biology. More importantly, we provide three extrapolation studies to demonstrate that PROSE-PDE can generalize physical features through the robust training of multiple operators and that the proposed model can extrapolate to predict PDE solutions whose models or data were unseen during the training. Furthermore, we show through systematic numerical experiments that the utilization of the symbolic modality in our model effectively resolves the well-posedness problems with training multiple operators and thus enhances our model's predictive capabilities.
We give a constructive characterization of matrices satisfying the reverse-order law for the Moore--Penrose pseudoinverse. In particular, for a given matrix $A$ we construct another matrix $B$, of arbitrary compatible size and chosen rank, in terms of the right singular vectors of $A$, such that the reverse order law for $AB$ is satisfied. Moreover, we show that any matrix satisfying this law comes from a similar construction. As a consequence, several equivalent conditions to $B^+ A^+$ being a pseudoinverse of $AB$ are given, for example $\mathcal{C}(A^*AB)=\mathcal{C}(BB^*A^*)$ or $B\left(AB\right)^+A$ being an orthogonal projection. In addition, we parameterize all possible SVD decompositions of a fixed matrix and give Greville-like equivalent conditions for $B^+A^+$ being a $\{1,2\}$-,$\{1,2,3\}$- and $\{1,2,4\}$-inverse of $AB$, with a geometric insight in terms of the principal angles between $\mathcal{C}(A^*)$ and $\mathcal{C}(B)$.
We introduce a simple and stable computational method for ill-posed partial differential equation (PDE) problems. The method is based on Schrödingerization, introduced in [S. Jin, N. Liu and Y. Yu, Phys. Rev. A, 108 (2023), 032603], which maps all linear PDEs into Schrödinger-type equations in one higher dimension, for quantum simulations of these PDEs. Although the original problem is ill-posed, the Schrödingerized equations are Hamiltonian systems and time-reversible, allowing stable computation both forward and backward in time. The original variable can be recovered by data from suitably chosen domain in the extended dimension. We will use the backward heat equation and the linear convection equation with imaginary wave speed as examples. Error analysis of these algorithms are conducted and verified numerically. The methods are applicable to both classical and quantum computers, and we also lay out quantum algorithms for these methods. Moreover, we introduce a smooth initialization for the Schrödingerized equation which will lead to essentially spectral accuracy for the approximation in the extended space, if a spectral method is used. Consequently, the extra qubits needed due to the extra dimension, if a qubit based quantum algorithm is used, for both well-posed and ill-posed problems, becomes almost $\log\log {1/\varepsilon}$ where $\varepsilon$ is the desired precision. This optimizes the complexity of the Schrödingerization based quantum algorithms for any non-unitary dynamical system introduced in [S. Jin, N. Liu and Y. Yu, Phys. Rev. A, 108 (2023), 032603].
Many interesting functions arising in applications map into Riemannian manifolds. We present an algorithm, using the manifold exponential and logarithm, for approximating such functions. Our approach extends approximation techniques for functions into linear spaces so that we can upper bound the forward error in terms of a lower bound on the manifold's sectional curvature. Furthermore, when the sectional curvature of a manifold is nonnegative, such as for compact Lie groups, the error is guaranteed to not be worse than in the linear case. We implement the algorithm in a Julia package and apply it to two example problems from Krylov subspaces and dynamic low-rank approximation, respectively. For these examples, the maps are confirmed to be well approximated by our algorithm.
In the present work we studied a subfield of Applied Mathematics called Riemannian Optimization. The main goal of this subfield is to generalize algorithms, theorems and tools from Mathematical Optimization to the case in which the optimization problem is defined on a Riemannian manifold. As a case study, we implemented some of the main algorithms described in the literature (Gradient Descent, Newton-Raphson and Conjugate Gradient) to solve an optimization problem known as Hartree-Fock. This method is extremely important in the field of Computational Quantum Chemistry and it is a good case study because it is a problem somewhat hard to solve and, as a consequence of this, it requires many tools from Riemannian Optimization. Besides, it is also a good example to see how these algorithms perform in practice.
In this paper, we investigate the reconstruction error, $N_\e^{\text{rec}}(x)$, when a linear, filtered back-projection (FBP) algorithm is applied to noisy, discrete Radon transform data with sampling step size $\epsilon$ in two-dimensions. Specifically, we analyze $N_\e^{\text{rec}}(x)$ for $x$ in small, $O(\e)$-sized neighborhoods around a generic fixed point, $x_0$, in the plane, where the measurement noise values, $\eta_{k,j}$ (i.e., the errors in the sinogram space), are random variables. The latter are independent, but not necessarily identically distributed. We show, under suitable assumptions on the first three moments of the $\eta_{k,j}$, that the following limit exists: $N^{\text{rec}}(\chx;x_0) = \lim_{\e\to0}N_\e^{\text{rec}}(x_0+\e\chx)$, for $\check x$ in a bounded domain. Here, $N_\e^{\text{rec}}$ and $ N^{\text{rec}}$ are viewed as continuous random variables, and the limit is understood in the sense of distributions. Once the limit is established, we prove that $N^{\text{rec}}$ is a zero mean Gaussian random field and compute explicitly its covariance. In addition, we validate our theory using numerical simulations and pseudo random noise.
Explicit, momentum-based dynamics for optimizing functions defined on Lie groups was recently constructed, based on techniques such as variational optimization and left trivialization. We appropriately add tractable noise to the optimization dynamics to turn it into a sampling dynamics, leveraging the advantageous feature that the trivialized momentum variable is Euclidean despite that the potential function lives on a manifold. We then propose a Lie-group MCMC sampler, by delicately discretizing the resulting kinetic-Langevin-type sampling dynamics. The Lie group structure is exactly preserved by this discretization. Exponential convergence with explicit convergence rate for both the continuous dynamics and the discrete sampler are then proved under $W_2$ distance. Only compactness of the Lie group and geodesically $L$-smoothness of the potential function are needed. To the best of our knowledge, this is the first convergence result for kinetic Langevin on curved spaces, and also the first quantitative result that requires no convexity or, at least not explicitly, any common relaxation such as isoperimetry.
The $k$-principal component analysis ($k$-PCA) problem is a fundamental algorithmic primitive that is widely-used in data analysis and dimensionality reduction applications. In statistical settings, the goal of $k$-PCA is to identify a top eigenspace of the covariance matrix of a distribution, which we only have black-box access to via samples. Motivated by these settings, we analyze black-box deflation methods as a framework for designing $k$-PCA algorithms, where we model access to the unknown target matrix via a black-box $1$-PCA oracle which returns an approximate top eigenvector, under two popular notions of approximation. Despite being arguably the most natural reduction-based approach to $k$-PCA algorithm design, such black-box methods, which recursively call a $1$-PCA oracle $k$ times, were previously poorly-understood. Our main contribution is significantly sharper bounds on the approximation parameter degradation of deflation methods for $k$-PCA. For a quadratic form notion of approximation we term ePCA (energy PCA), we show deflation methods suffer no parameter loss. For an alternative well-studied approximation notion we term cPCA (correlation PCA), we tightly characterize the parameter regimes where deflation methods are feasible. Moreover, we show that in all feasible regimes, $k$-cPCA deflation algorithms suffer no asymptotic parameter loss for any constant $k$. We apply our framework to obtain state-of-the-art $k$-PCA algorithms robust to dataset contamination, improving prior work in sample complexity by a $\mathsf{poly}(k)$ factor.
We reformulate the Lanczos tau method for the discretization of time-delay systems in terms of a pencil of operators, allowing for new insights into this approach. As a first main result, we show that, for the choice of a shifted Legendre basis, this method is equivalent to Padé approximation in the frequency domain. We illustrate that Lanczos tau methods straightforwardly give rise to sparse, self nesting discretizations. Equivalence is also demonstrated with pseudospectral collocation, where the non-zero collocation points are chosen as the zeroes of orthogonal polynomials. The importance of such a choice manifests itself in the approximation of the $H^2$-norm, where, under mild conditions, super-geometric convergence is observed and, for a special case, super convergence is proved; both significantly faster than the algebraic convergence reported in previous work.
We study symmetric tensor decompositions, i.e., decompositions of the form $T = \sum_{i=1}^r u_i^{\otimes 3}$ where $T$ is a symmetric tensor of order 3 and $u_i \in \mathbb{C}^n$.In order to obtain efficient decomposition algorithms, it is necessary to require additional properties from $u_i$. In this paper we assume that the $u_i$ are linearly independent. This implies $r \leq n$,that is, the decomposition of T is undercomplete. We give a randomized algorithm for the following problem in the exact arithmetic model of computation: Let $T$ be an order-3 symmetric tensor that has an undercomplete decomposition.Then given some $T'$ close to $T$, an accuracy parameter $\varepsilon$, and an upper bound B on the condition number of the tensor, output vectors $u'_i$ such that $||u_i - u'_i|| \leq \varepsilon$ (up to permutation and multiplication by cube roots of unity) with high probability. The main novel features of our algorithm are: 1) We provide the first algorithm for this problem that runs in linear time in the size of the input tensor. More specifically, it requires $O(n^3)$ arithmetic operations for all accuracy parameters $\varepsilon =$ 1/poly(n) and B = poly(n). 2) Our algorithm is robust, that is, it can handle inverse-quasi-polynomial noise (in $n$,B,$\frac{1}{\varepsilon}$) in the input tensor. 3) We present a smoothed analysis of the condition number of the tensor decomposition problem. This guarantees that the condition number is low with high probability and further shows that our algorithm runs in linear time, except for some rare badly conditioned inputs. Our main algorithm is a reduction to the complete case ($r=n$) treated in our previous work [Koiran,Saha,CIAC 2023]. For efficiency reasons we cannot use this algorithm as a blackbox. Instead, we show that it can be run on an implicitly represented tensor obtained from the input tensor by a change of basis.
In this work, we analyze a sublinear-time algorithm for selecting a few rows and columns of a matrix for low-rank approximation purposes. The algorithm is based on an initial uniformly random selection of rows and columns, followed by a refinement of this choice using a strong rank-revealing QR factorization. We prove bounds on the error of the corresponding low-rank approximation (more precisely, the CUR approximation error) when the matrix is a perturbation of a low-rank matrix that can be factorized into the product of matrices with suitable incoherence and/or sparsity assumptions.
We reiterate the contribution made by Harrow, Hassidim, and Llyod to the quantum matrix equation solver with the emphasis on the algorithm description and the error analysis derivation details. Moreover, the behavior of the amplitudes of the phase register on the completion of the Quantum Phase Estimation is studied. This study is beneficial for the comprehension of the choice of the phase register size and its interrelation with the Hamiltonian simulation duration in the algorithm setup phase.
Chebyshev varieties are algebraic varieties parametrized by Chebyshev polynomials or their multivariate generalizations. We determine the dimension, degree, singular locus and defining equations of these varieties. We explain how they play the role of toric varieties in sparse polynomial root finding, when monomials are replaced by Chebyshev polynomials. We present numerical root finding algorithms that exploit our results.
This report is concerned with the efficiency of numerical methods for simulating quantum spin systems, with the aim to implement an improved method for simulation of a time-dependent Hamiltonian that displays chirped pulses at a high frequency. Working in the density matrix formulation of quantum systems, we study evolution under the Liouville-von Neumann equation, presenting analysis of and benchmarking current numerical methods. The accuracy of existing techniques is assessed in the presence of chirped pulses. We also discuss the Magnus expansion and detail how a truncation of it is used to solve differential equations. The results of this work are implemented in the Python package MagPy to provide a better error-to-cost ratio than current approaches allow for time-dependent Hamiltonians.
One of the pivotal tasks in scientific machine learning is to represent underlying dynamical systems from time series data. Many methods for such dynamics learning explicitly require the derivatives of state data, which are not directly available and can be approximated conventionally by finite differences. However, the discrete approximations of time derivatives may result in poor estimations when state data are scarce and/or corrupted by noise, thus compromising the predictiveness of the learned dynamical models. To overcome this technical hurdle, we propose a new method that learns nonlinear dynamics through a Bayesian inference of characterizing model parameters. This method leverages a Gaussian process representation of states, and constructs a likelihood function using the correlation between state data and their derivatives, yet prevents explicit evaluations of time derivatives. Through a Bayesian scheme, a probabilistic estimate of the model parameters is given by the posterior distribution, and thus a quantification is facilitated for uncertainties from noisy state data and the learning process. Specifically, we will discuss the applicability of the proposed method to several typical scenarios for dynamical systems: identification and estimation with an affine parametrization, nonlinear parametric approximation without prior knowledge, and general parameter estimation for a given dynamical system.
In this chapter we provide a theoretically founded investigation of state-of-the-art learning approaches for inverse problems from the point of view of spectral reconstruction operators. We give an extended definition of regularization methods and their convergence in terms of the underlying data distributions, which paves the way for future theoretical studies. Based on a simple spectral learning model previously introduced for supervised learning, we investigate some key properties of different learning paradigms for inverse problems, which can be formulated independently of specific architectures. In particular we investigate the regularization properties, bias, and critical dependence on training data distributions. Moreover, our framework allows to highlight and compare the specific behavior of the different paradigms in the infinite-dimensional limit.
We give a stochastic optimization algorithm that solves a dense $n\times n$ real-valued linear system $Ax=b$, returning $\tilde x$ such that $\|A\tilde x-b\|\leq \epsilon\|b\|$ in time: $$\tilde O((n^2+nk^\omega-1)\log1/\epsilon),$$ where $k$ is the number of singular values of $A$ larger than $O(1)$ times its smallest positive singular value, $\omega < 2.372$ is the matrix multiplication exponent, and $\tilde O$ hides a poly-logarithmic in $n$ factor. When $k=O(n^{1-\theta})$ (namely, $A$ has a flat-tailed spectrum, e.g., due to noisy data or regularization), this improves on both the cost of solving the system directly, as well as on the cost of preconditioning an iterative method such as conjugate gradient. In particular, our algorithm has an $\tilde O(n^2)$ runtime when $k=O(n^{0.729})$. We further adapt this result to sparse positive semidefinite matrices and least squares regression. Our main algorithm can be viewed as a randomized block coordinate descent method, where the key challenge is simultaneously ensuring good convergence and fast per-iteration time. In our analysis, we use theory of majorization for elementary symmetric polynomials to establish a sharp convergence guarantee when coordinate blocks are sampled using a determinantal point process. We then use a Markov chain coupling argument to show that similar convergence can be attained with a cheaper sampling scheme, and accelerate the block coordinate descent update via matrix sketching.
This work is concerned with solving high-dimensional Fokker-Planck equations with the novel perspective that solving the PDE can be reduced to independent instances of density estimation tasks based on the trajectories sampled from its associated particle dynamics. With this approach, one sidesteps error accumulation occurring from integrating the PDE dynamics on a parameterized function class. This approach significantly simplifies deployment, as one is free of the challenges of implementing loss terms based on the differential equation. In particular, we introduce a novel class of high-dimensional functions called the functional hierarchical tensor (FHT). The FHT ansatz leverages a hierarchical low-rank structure, offering the advantage of linearly scalable runtime and memory complexity relative to the dimension count. We introduce a sketching-based technique that performs density estimation over particles simulated from the particle dynamics associated with the equation, thereby obtaining a representation of the Fokker-Planck solution in terms of our ansatz. We apply the proposed approach successfully to three challenging time-dependent Ginzburg-Landau models with hundreds of variables.
Low-rank approximation of a matrix function, $f(A)$, is an important task in computational mathematics. Most methods require direct access to $f(A)$, which is often considerably more expensive than accessing $A$. Persson and Kressner (SIMAX 2023) avoid this issue for symmetric positive semidefinite matrices by proposing funNyström, which first constructs a Nyström approximation to $A$ using subspace iteration, and then uses the approximation to directly obtain a low-rank approximation for $f(A)$. They prove that the method yields a near-optimal approximation whenever $f$ is a continuous operator monotone function with $f(0) = 0$. We significantly generalize the results of Persson and Kressner beyond subspace iteration. We show that if $\widehat{A}$ is a near-optimal low-rank Nyström approximation to $A$ then $f(\widehat{A})$ is a near-optimal low-rank approximation to $f(A)$, independently of how $\widehat{A}$ is computed. Further, we show sufficient conditions for a basis $Q$ to produce a near-optimal Nyström approximation $\widehat{A} = AQ(Q^T AQ)^{\dagger} Q^T A$. We use these results to establish that many common low-rank approximation methods produce near-optimal Nyström approximations to $A$ and therefore to $f(A)$.
We propose a new algorithm for efficiently solving the damped Fisher matrix in large-scale scenarios where the number of parameters significantly exceeds the number of available samples. This problem is fundamental for natural gradient descent and stochastic reconfiguration. Our algorithm is based on Cholesky decomposition and is generally applicable. Benchmark results show that the algorithm is significantly faster than existing methods.
Nonlinear boolean equation systems play an important role in a wide range of applications. Grover's algorithm is one of the best-known quantum search algorithms in solving the nonlinear boolean equation system on quantum computers. In this paper, we propose three novel techniques to improve the efficiency under Grover's algorithm framework. A W-cycle circuit construction introduces a recursive idea to increase the solvable number of boolean equations given a fixed number of qubits. Then, a greedy compression technique is proposed to reduce the oracle circuit depth. Finally, a randomized Grover's algorithm randomly chooses a subset of equations to form a random oracle every iteration, which further reduces the circuit depth and the number of ancilla qubits. Numerical results on boolean quadratic equations demonstrate the efficiency of the proposed techniques.
For a reaction-dominated diffusion problem we study a primal and a dual hybrid finite element method where weak continuity conditions are enforced by Lagrange multipliers. Uniform robustness of the discrete methods is achieved by enriching the local discretization spaces with modified face bubble functions which decay exponentially in the interior of an element depending on the ratio of the singular perturbation parameter and the local mesh-size. A posteriori error estimators are derived using Fortin operators. They are robust with respect to the singular perturbation parameter. Numerical experiments are presented that show that oscillations, if present, are significantly smaller then those observed in common finite element methods.
We propose nodal auxiliary space preconditioners for facet and edge virtual elements of lowest order by deriving discrete regular decompositions on polytopal grids and generalizing the Hiptmair-Xu preconditioner to the virtual element framework. The preconditioner consists of solving a sequence of elliptic problems on the nodal virtual element space, combined with appropriate smoother steps. Under assumed regularity of the mesh, the preconditioned system is proven to have bounded spectral condition number independent of the mesh size and this is verified by numerical experiments on a sequence of polygonal meshes. Moreover, we observe numerically that the preconditioner is robust on meshes containing elements with high aspect ratios.
This work presents a new algorithm to compute the matrix exponential within a given tolerance. Combined with the scaling and squaring procedure, the algorithm incorporates Taylor, partitioned and classical Padé methods shown to be superior in performance to the approximants used in state-of-the-art software. The algorithm computes matrix--matrix products and also matrix inverses, but it can be implemented to avoid the computation of inverses, making it convenient for some problems. If the matrix A belongs to a Lie algebra, then exp(A) belongs to its associated Lie group, being a property which is preserved by diagonal Padé approximants, and the algorithm has another option to use only these. Numerical experiments show the superior performance with respect to state-of-the-art implementations.
The conjugate function method is an algorithm for numerical computation of conformal mappings for simply and multiply connected domains. In this paper, the conjugate function method is extended to cover conformal mappings between Riemannian surfaces. The main challenge addressed here is the connection between Laplace--Beltrami equations on surfaces and the computation of the conformal modulus of a quadrilateral. We consider mappings of simply, doubly, and multiply connected domains. The numerical computation is based on an $hp$-adaptive finite element method. The key advantage of our approach is that it allows highly accurate computations of mappings on surfaces, including domains of complex boundary geometry involving strong singularities and cusps. The efficacy of the proposed method is illustrated via an extensive set of numerical experiments including error estimates.
In this work we study the stability, convergence, and pressure-robustness of discretization methods for incompressible flows with hybrid velocity and pressure. Specifically, focusing on the Stokes problem, we identify a set of assumptions that yield inf-sup stability as well as error estimates which distinguish the velocity- and pressure-related contributions to the error. We additionally identify the key properties under which the pressure-related contributions vanish in the estimate of the velocity, thus leading to pressure-robustness. Several examples of existing and new schemes that fit into the framework are provided, and extensive numerical validation of the theoretical properties is provided.
This paper proposes an effective treatment of hyperparameters in the Bayesian inference of a scalar field from indirect observations. Obtaining the joint posterior distribution of the field and its hyperparameters is challenging. The infinite dimensionality of the field requires a finite parametrization that usually involves hyperparameters to reflect the limited prior knowledge. In the present work, we consider a Karhunen-Loève (KL) decomposition for the random field and hyperparameters to account for the lack of prior knowledge of its autocovariance function. The hyperparameters must be inferred. To efficiently sample jointly the KL coordinates of the field and the autocovariance hyperparameters, we introduce a change of measure to reformulate the joint posterior distribution into a hierarchical Bayesian form. The likelihood depends only on the field's coordinates in a fixed KL basis, with a prior conditioned on the hyperparameters. We exploit this structure to derive an efficient Markov Chain Monte Carlo (MCMC) sampling scheme based on an adapted Metropolis-Hasting algorithm. We rely on surrogate models (Polynomial Chaos expansions) of the forward model predictions to further accelerate the MCMC sampling. A first application to a transient diffusion problem shows that our method is consistent with other approaches based on a change of coordinates (Sraj et al., 2016). A second application to a seismic traveltime tomography highlights the importance of inferring the hyperparameters.
In a Jacobi--Davidson (JD) type method for singular value decomposition (SVD) problems, called JDSVD, a large symmetric and generally indefinite correction equation is approximately solved iteratively at each outer iteration, which constitutes the inner iterations and dominates the overall efficiency of JDSVD. In this paper, a convergence analysis is made on the minimal residual (MINRES) method for the correction equation. Motivated by the results obtained, a preconditioned correction equation is derived that extracts useful information from current searching subspaces to construct effective preconditioners for the correction equation and is proved to retain the same convergence of outer iterations of JDSVD. The resulting method is called inner preconditioned JDSVD (IPJDSVD) method. Convergence results show that MINRES for the preconditioned correction equation can converge much faster when there is a cluster of singular values closest to a given target, so that IPJDSVD is more efficient than JDSVD. A new thick-restart IPJDSVD algorithm with deflation and purgation is proposed that simultaneously accelerates the outer and inner convergence of the standard thick-restart JDSVD and computes several singular triplets of a large matrix. Numerical experiments justify the theory and illustrate the considerable superiority of IPJDSVD to JDSVD.
In this paper, we first propose a coupled numerical model of unsaturated flow in soils and plant root water uptake. The Richards equation and different formulations are used in the developed numerical model to describe infiltration in root zone and to investigate the impact of the plant root on the distribution of soil moisture. The Kirchhoff transformed Richards equation is used and the Gardner model is considered for capillary pressure. In our approach, we employ a meshless method based on localized radial basis functions (LRBF) to solve the resulting system of equations. The LRBF approach is an accurate and computationally efficient method that does not require mesh generation and is flexible in addressing high-dimensional problems with complex geometries. Furthermore, this method leads to a sparse matrix system, which avoids ill-conditioning issues. We implement the coupled numerical model of infiltration and plant root water uptake for one, two, and three-dimensional soils. Numerical experiments are performed using nontrivial analytical solutions and available experimental data to validate the coupled numerical model. The numerical results demonstrate the performance and ability of the proposed numerical method to predict soil moisture dynamics in root zone.
The demagnetization field in micromagnetism is given as the gradient of a potential which solves a partial differential equation (PDE) posed in R^d. In its most general form, this PDE is supplied with continuity condition on the boundary of the magnetic domain and the equation includes a discontinuity in the gradient of the potential over the boundary. Typical numerical algorithms to solve this problem relies on the representation of the potential via the Green's function, where a volume and a boundary integral terms need to be accurately approximated. From a computational point of view, the volume integral dominates the computational cost and can be difficult to approximate due to the singularities of the Green's function. In this article, we propose a hybrid model, where the overall potential can be approximated by solving two uncoupled PDEs posed in bounded domains, whereby the boundary conditions of one of the PDEs is obtained by a low cost boundary integral. Moreover, we provide a convergence analysis of the method under two separate theoretical settings; periodic magnetisation, and high-frequency magnetisation. Numerical examples are given to verify the convergence rates.
This paper is concerned with model order reduction of parametric Partial Differential Equations (PDEs) using tree-based library approximations. Classical approaches are formulated for PDEs on Hilbert spaces and involve one single linear space to approximate the set of PDE solutions. Here, we develop reduced models relying on a collection of linear or nonlinear approximation spaces called a library, and which can also be formulated on general metric spaces. To build the spaces of the library, we rely on greedy algorithms involving different splitting strategies which lead to a hierarchical tree-based representation. We illustrate through numerical examples that the proposed strategies have a much wider range of applicability in terms of the parametric PDEs that can successfully be addressed. While the classical approach is very efficient for elliptic problems with strong coercivity, we show that the tree-based library approaches can deal with diffusion problems with weak coercivity, convection-diffusion problems, and with transport-dominated PDEs posed on general metric spaces such as the $L^2$-Wasserstein space.
In this work we provide theoretical estimates for the ranks of the power functions $f(k) = k^{-\alpha}$, $\alpha>1$ in the quantized tensor train (QTT) format for $k = 1, 2, 3, \ldots, 2^{d}$. Such functions and their several generalizations (e.~g. $f(k) = k^{-\alpha} \cdot e^{-\lambda k}, \lambda > 0$) play an important role in studies of the asymptotic solutions of the aggregation-fragmentation kinetic equations. In order to support the constructed theory we verify the values of QTT-ranks of these functions in practice with the use of the TTSVD procedure and show an agreement between the numerical and analytical results.
Parallelisation in Bayesian optimisation is a common strategy but faces several challenges: the need for flexibility in acquisition functions and kernel choices, flexibility dealing with discrete and continuous variables simultaneously, model misspecification, and lastly fast massive parallelisation. To address these challenges, we introduce a versatile and modular framework for batch Bayesian optimisation via probabilistic lifting with kernel quadrature, called SOBER, which we present as a Python library based on GPyTorch/BoTorch. Our framework offers the following unique benefits: (1) Versatility in downstream tasks under a unified approach. (2) A gradient-free sampler, which does not require the gradient of acquisition functions, offering domain-agnostic sampling (e.g., discrete and mixed variables, non-Euclidean space). (3) Flexibility in domain prior distribution. (4) Adaptive batch size (autonomous determination of the optimal batch size). (5) Robustness against a misspecified reproducing kernel Hilbert space. (6) Natural stopping criterion.
In the context of the optimization of rotating electric machines, many different objective functions are of interest and considering this during the optimization is of crucial importance. While evolutionary algorithms can provide a Pareto front straightforwardly and are widely used in this context, derivative-based optimization algorithms can be computationally more efficient. In this case, a Pareto front can be obtained by performing several optimization runs with different weights. In this work, we focus on a free-form shape optimization approach allowing for arbitrary motor geometries. In particular, we propose a way to efficiently obtain Pareto-optimal points by moving along to the Pareto front exploiting a homotopy method based on second order shape derivatives.
Braneshop, and
the US National Science Foundation.
SciRate