allowed lowdin with SVD and to return the transformation matrix

Signed-off-by: Nick Papior <[email protected]>

CHANGELOG.md

@@ -8,6 +8,9 @@ we hit release version 1.0.0.
 
 ## [0.15.1] - YYYY-MM-DD
 
+### Added
+- enabled `lowdin` to return the Lowdin transformation matrix, and also
+  allow it to be calculated using SVD
 
 
 ## [0.15.0] - 2024-08-13

src/sisl/linalg/special.py

@@ -3,9 +3,11 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 from __future__ import annotations
 
+from typing import Literal
+
 import numpy as np
 
-from .base import eigh
+from .base import eigh, svd
 
 __all__ = ["signsqrt", "sqrth", "invsqrth", "lowdin"]
 
@@ -53,19 +55,48 @@ def invsqrth(a, overwrite_a=False):
     return (ev * eig) @ ev.conj().T
 
 
-def lowdin(a, b, overwrite_a=False):
-    r"""Convert the matrix `b` in the basis `a` into an orthogonal basis using the Lowdin transformation
+def lowdin(
+    a,
+    b=None,
+    overwrite_a: bool = False,
+    driver: Literal["eigh", "gesdd", "gesvd"] = "eigh",
+):
+    r"""Calculate the Lowdin transformation matrix, optionally convert the matrix `b` into the orthogonal basis
+
+    Convert the matrix `b` in the basis `a` into an orthogonal basis using the Lowdin transformation
 
     .. math::
 
        \mathbf B' = \mathbf A^{-1/2} \mathbf B \mathbf A^{-1/2}
 
+    Note, this method assumes `a` to be Hermitian.
+
     Parameters
     ----------
     a : array_like
        basis matrix to convert to the Lowdin basis
     b : array_like
-       matrix to convert
+       matrix to convert, if not provided as an argument, :math:`\mathbf A^{-1/2}` is
+       returned.
+    overwrite_a :
+        specificy whether `a` can be altered in the call.
+    driver :
+        which driver to use for calculating the Lowdin transformed `a` matrix.
+        When `a` is a matrix with rank ``len(a)`` with algebraic multiplicity
+        :math:`\mu_A(1)` equal to the rank, then the SVD method can be used
+        to construct the Lowdin transformation.
     """
-    a12 = invsqrth(a, overwrite_a=overwrite_a)
+    if driver == "eigh":
+        a12 = invsqrth(a, overwrite_a=overwrite_a)
+    elif driver.startswith("ges"):
+        U, _, Vh = svd(
+            a, compute_uv=True, overwrite_a=overwrite_a, lapack_driver=driver
+        )
+        a12 = U @ Vh
+        del U, _, Vh
+    else:
+        raise ValueError(f"lowdin: got unknown driver argument '{driver}'")
+
+    if b is None:
+        return a12
     return a12 @ b @ a12

src/sisl/linalg/tests/test_special.py

@@ -7,7 +7,7 @@
 import pytest
 import scipy.linalg as sl
 
-from sisl.linalg import invsqrth, signsqrt, sqrth
+from sisl.linalg import invsqrth, lowdin, signsqrt, sqrth
 
 pytestmark = [pytest.mark.linalg]
 
@@ -61,3 +61,32 @@ def test_invsqrth_offset():
     np.divide(1, eig, where=(eig != 0), out=eig)
     sa2 = (ev * eig) @ ev.conj().T
     assert np.allclose(sa1, sa2)
+
+
+@pytest.mark.parametrize("driver", ["eigh", "gesdd", "gesvd"])
+def test_lowdin(driver):
+    # offsetting eigenvalues only works if the matrix is
+    # positive semi-definite
+    np.random.seed(1285947159)
+
+    b = np.random.rand(10, 10)
+    b = b + b.T
+    np.fill_diagonal(b, b.diagonal() + 4)
+
+    # Recreate an overlap matrix with all eigenvalues being 1.
+    _, ev = sl.eigh(b)
+    b = ev.T @ ev
+
+    # Create "H"
+    a = np.random.rand(10, 10)
+    a = a + a.T
+
+    eig, ev = sl.eigh(a, b)
+
+    # Calculate eigenvalues for a, b and then
+    # by the Lowdin transformed a
+    aL = lowdin(b, a, driver=driver)
+    eigL, evL = sl.eigh(aL)
+
+    assert np.allclose(eig, eigL)
+    assert not np.allclose(ev, evL)