Overlap matrix and its inverse matrix

Question

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.

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.