# Overlap matrix and its inverse matrix

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.

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

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.