Now, we consider a non-orthonormal basis: $$\mathcal{S}_N=\{|\alpha\rangle,a^\dagger|\alpha\rangle,a^{\dagger 2}|\alpha\rangle,\ldots,a^{\dagger N}|\alpha\rangle\},$$ where $|\alpha\rangle$ is the coherent state and $a$ is the annihilation operator of bosonic mode.

Then, we assume $|\phi_m\rangle=a^{\dagger m}|\alpha\rangle\in\mathcal{S}_N$, and define the overlap matrix $S$ with the matrix elements: $$S_{n,m}=\langle\phi_n|\phi_m\rangle, \text{where}\quad n,m\in[0,N].$$

In general, the problem of calculating the overlap matrix $S$ is simple, but its inverse matrix is not easy.

Finally, even if we can't obtain the analytical expression about the inverse matrix, I also want to obtain it efficiently in Matlab or Python. I want to explain why we can't correctly obtain its inverse matrix in the program. Firstly, when $|\alpha|\gg1$, the matrix $S$ becomes an ill-matrix, so we can't correctly obtain its inverse matrix, i.e., $S^{-1}S\neq I$, or the error is huge when we use the inverse matrix to do matrix multiplication. Secondly, we also choose a symbolic language to solve it, but the price greatly increases the computation time of matrix multiplication and addition. Finally, I especially want to find a suitable algorithm to solve it, and the analytical expression is secondary because sometimes I need to replace $|\alpha\rangle$ with $[|\alpha\rangle\pm|-\alpha\rangle]$ in my basis.

  • 1
    $\begingroup$ It is commonly adviced to not compute inverse of matrices. Most of the times you can solve a linear system instead. $\endgroup$
    – nicoguaro
    Commented Dec 4, 2023 at 2:25
  • 2
    $\begingroup$ What are typical values for $N$ and for the length of vectors involved? Is $a$ sparse? Does the dagger stand for transposition? $\endgroup$ Commented Dec 4, 2023 at 7:31
  • $\begingroup$ @FedericoPoloni 1<N<30, the length of vectors is N+1, a is sparse, and the dagger stand for conjugate and transpose. $\endgroup$
    – Young Q
    Commented Dec 5, 2023 at 7:20

1 Answer 1


There are several improvements you can make to the computation of $S_{n,m}^{-1}$ to make it more stable.

Avoid the explicit inverse

As the comments say, this is the first thing to look into. You should aim to compute a decomposition of the matrix as a product of "easy-to-invert" factors (triangular, orthogonal, permutations). Then you can solve any linear system with your matrix by computing the action of the inverse of those factors on a vector, one by one: for instance, if you find a factorization $A=QR$, with $Q$ unitary and $R$ upper triangular, you can solve the linear system $b = Ax= QRx$ by computing first $y = Q^\dagger b$ and then solving $Rx=y$ by back-substitution.

This is the standard way to work with linear systems and matrix inverses in practice. Look up the LU and QR factorization (not the QR iteration / algorithm, which is another thing) to get more information on this.

Avoid forming $A^\dagger A$

Define $$ A=[|\alpha\rangle,a^\dagger|\alpha\rangle,a^{\dagger 2}|\alpha\rangle,\ldots,a^{\dagger N}|\alpha\rangle] \in \mathbb{C}^{(N+1)\times (N+1)}. $$ You have $S_{n,m} = A^\dagger A$, so $S_{n,m}^{-1} = A^{-1} (A^\dagger)^{-1}$. The second improvement is looking for a factorization of $A$ rather than one of $S_{n,m}$. The matrix $S_{n,m}$ is more ill-conditioned than $A$, since $\kappa(S_{n,m}) = \kappa(A)^2$, so solving a linear system with $S_{n,m}$ gives a larger error than solving one with $A$ and one with $A^\dagger $, in sequence.

Avoid forming $A$ using Arnoldi

The third improvement is not constructing $A$. The matrix $A$ may be very ill-conditioned by itself, since its columns tend to the leading eigenvector of $A$ (the one with the largest associated eigenvalue, in modulus). Look up the power method to find out more about why.

If you first compute $A$ and then factorize it, your methods cannot have a better relative error than $\kappa(A)$. Fortunately, there is a way to compute a certain factorization of the matrix $A$ in your problem without even forming it: the Arnoldi algorithm.

The Arnoldi algorithm (which you can run up to the size of the matrix if $N\approx 30$) gives you a factorization $A = QHQ^\dagger$, with $Q$ unitary and $H$ a Hessenberg matrix (i.e., $H_{ij}=0$ if $i>j+1$).

There are tricks to convert this factorization into a more standard QR factorization, but actually you can just use it to solve linear systems as it is, as solving linear systems with Hessenberg matrices costs $O(N^2)$. Essentially you have to compute a QR factorization of $H$ using Givens rotations. Unfortunately this last algorithm is a little more technical and not readily available in Lapack, for instance.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.