# How to do Lowdin symmetric orthonormalisation?

Löwdin symmetric orthonormalisation seems to be a common practice in quantum chemistry. I come from a different background though and have to understand it and possibly implement in a computer code. Here is the definition of the transformation;

$$|{\phi'}\rangle = S^{-1/2}|\phi\rangle,$$

where $S$ is the overlap matrix between the non-orthogonal basis,

$$S = \left< \phi|\phi \right> .$$

1. My first question is how to obtain $S^{-1/2}$ analytically?

2. Are there any parallel algorithms available to compute it numerically?

• I would recommend checking out Szabo and Ostlund Modern Quantum Chemistry. Its a very affordable book that discusses a lot of the framework of computational chemistry.
– Tyberius
Nov 8, 2017 at 17:34

You obtain the diagonalization of $S$ as defined $$S = U \cdot s \cdot U^{\dagger}$$ (which holds for the case of symmetric $S$, which it is) and apply the power needed $$S^{-\frac{1}{2}} = U \cdot s^{-\frac{1}{2}} \cdot U^{\dagger}$$

For more background, try the application section of Wikipedia's article on diagonalizable matrices. As far as understanding it goes, I'm afraid not much beyond doing some model calculations on 2D or 3D vectors can be done.

I do not have a strong background in parallel computing, but all operations involved are well-known and your linear algebra framework (BLAS, anyone?) will document how to run any of this in parallel and what storage may be needed for it. Note that, typically, diagonalization is not the most expensive step of a self-consistent field calculation (Hartree-Fock or DFT) anyway, so parallelization may not be critical for this step.

• Hi! I know that this was posted a while ago, but I would appreciate your help. How did you know that you diagonalize S as UsU^T? Aug 22, 2018 at 13:31
• @JunJang That is the basic property of a diagonalizable, symmetric matrix. I would not know what information to give you besides the linked wiki article or to point you to linear algebra textbooks. Aug 23, 2018 at 4:45
• I see. Will definitely look into wiki/textbooks. May I ask you one more question? If you were to explain the lowdin orthogonalization to a 5th grade, how would you explain it? I am working in finance and worked on a project using the lowdin orthogonalization and have to present it to management. I am using this to identify the underlying uncorrelated components of the factors and maintains the interpretations of the original factors.. How would you describe the orthogonal transformation and the associated diagonalization step in layman's terms? Thank you very much!!!! Aug 23, 2018 at 11:58
• Not sure that can be done easily, any decent explanation would have to resort to Taylor series. There are illustrations on how vectors are changed by orthogonalization, but I cannot find a good one. Aug 23, 2018 at 20:33
• would there be any way to explain the symmetric orthogonalization to someone without sufficient mathematics background? Aug 28, 2018 at 12:28

Your first question has already been answered, but in words: to find $f(A)$ for some matrix $A$, you diagonalize it to obtain the eigenvalues $a$ and eigenvectors $U$, apply $f$ to the diagonalized matrix (the eigenvalues), then back-transform $f(a)$ using the eigenvectors to the original non-diagonal basis.

For the second question, yes. Most revolve around calling the LAPACK routine dsyev (real double-precision symmetric eigen decomposition) or its variants, ssyev for real single-precision, {c,z}heev for complex single- and double-precision Hermitian decomposition. The vast majority of quantum chemistry uses entirely real coefficients and double-precision.

Assume I have already compute the overlap matrix $S$ by some method, usually by calling a program's integral engine. Here is a sample implementation in Python using NumPy:

print("Overlap matrix")
print(S)

lam_s, l_s = np.linalg.eigh(S)
lam_s = lam_s * np.eye(len(lam_s))
lam_sqrt_inv = np.sqrt(np.linalg.inv(lam_s))
symm_orthog = np.dot(l_s, np.dot(lam_sqrt_inv, l_s.T))

print("Symmetric orthogonalization matrix")
print(symm_orthog)


From the documentation of numpy.linalg.eigh:

The eigenvalues/eigenvectors are computed using LAPACK routines _syevd, _heevd

Here is a sample implementation in C++ using Armadillo:

S.print("Overlap matrix");

arma::vec lam_s_vec;
arma::mat l_s;
arma::eig_sym(lam_s_vec, l_s, S);
arma::mat lam_s_mat = arma::diagmat(lam_s_vec);
arma::mat lam_sqrt_inv = arma::sqrt(arma::inv(lam_s_mat));
arma::mat symm_orthog = l_s * lam_sqrt_inv * l_s.t();

symm_orthog.print("Symmetric orthogonalization matrix");


Both working examples can be found here; check the Makefile for how to run.

In the Armadillo source tree, the files include/armadillo_bits/{def,wrapper}_lapack.hpp contain more information about which LAPACK routines are called for which types.

Regarding parallelization, your BLAS + LAPACK implementation (MKL, OpenBLAS, ATLAS, ...) is most likely threaded and can be controlled by <something>_NUM_THREADS=4, where <something> might be OMP, MKL, OpenBLAS, or possibly something else, but you should check the documentation. This means that as long as your environment is set up properly, math library calls with NumPy or C++ template libraries like Armadillo or Eigen will run in parallel without explicit OpenMP annotations or MPI code. For distributed parallelization (MPI), there is ScaLAPACK, which shares a similar interface to regular LAPACK.

• Just to compare to a possible alternative, scipy.linalg.fractional_matrix_power(S, -0.5) gives the same result, but is much slower (~15x slower for a 24 x 24 overlap matrix) than your approach :) Aug 29, 2019 at 2:38
• You do not have to compute the inverse. The matrix is already diagonalized. Just compute the power of the vector with the eigenvalues and then construct the matrix with the vector as the main diagonal. That is, np.diag(lam_s**(-0.5)). Jul 18, 2020 at 11:08
• Also, symm_ortho = np.einsum("ik,kj,lj->il",lam_s_vec, lam_sqrt_inv, lam_s_vec) could be faster. Jul 18, 2020 at 11:13