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.
