LU instead of L

.gitignore

@@ -1,3 +1,4 @@
+poetry.lock
 .pytest_cache/
 **/__pycache__/
 .vscode

.pre-commit-config.yaml

@@ -1,19 +1,15 @@
 repos:
     - repo: https://github.com/ambv/black
-      rev: 21.5b1
+      rev: 22.6.0
       hooks:
           - id: black
             language_version: python3
     - repo: https://github.com/timothycrosley/isort
-      rev: 5.8.0
+      rev: 5.10.1
       hooks:
           - id: isort
-    - repo: https://gitlab.com/pycqa/flake8
-      rev: 3.9.2
-      hooks:
-          - id: flake8
     - repo: https://github.com/pre-commit/pre-commit-hooks
-      rev: v4.0.1
+      rev: v4.3.0
       hooks:
           - id: check-added-large-files
           - id: trailing-whitespace

glimix_core/__init__.py

@@ -6,7 +6,7 @@
 from . import cov, example, ggp, glmm, gp, lik, link, lmm, mean, random
 from ._testit import test
 
-__version__ = "3.1.12"
+__version__ = "3.1.13"
 
 __all__ = [
     "__version__",

glimix_core/_ep/ep.py

@@ -122,14 +122,14 @@ def cav(self):
         return self._cav
 
     def lml(self):
-        from numpy_sugar.linalg import cho_solve
+        from numpy_sugar.linalg import lu_slogdet, lu_solve
 
         if self._cache["lml"] is not None:
             return self._cache["lml"]
 
         self._update()
 
-        L = self._posterior.L()
+        LU = self._posterior.LU()
         Q, S = self._posterior.cov["QS"]
         Q = Q[0]
         ttau = self._site.tau
@@ -141,13 +141,13 @@ def lml(self):
         TS = ttau + ctau
 
         lml = [
-            -log(L.diagonal()).sum(),
+            -lu_slogdet(LU)[1],
             -0.5 * sum(log(S)),
             # lml += 0.5 * sum(log(ttau)),
-            +0.5 * dot(teta, dot(Q, cho_solve(L, dot(Q.T, teta)))),
+            +0.5 * dot(teta, dot(Q, lu_solve(LU, dot(Q.T, teta)))),
             -0.5 * dot(teta, teta / TS),
             +dot(m, teta) - 0.5 * dot(m, ttau * m),
-            -0.5 * dot(m * ttau, dot(Q, cho_solve(L, dot(Q.T, 2 * teta - ttau * m)))),
+            -0.5 * dot(m * ttau, dot(Q, lu_solve(LU, dot(Q.T, 2 * teta - ttau * m)))),
             +sum(self._moments["log_zeroth"]),
             +0.5 * sum(log(TS)),
             # lml -= 0.5 * sum(log(ttau)),
@@ -163,45 +163,45 @@ def lml(self):
         return lml
 
     def lml_derivative_over_cov(self, dQS):
-        from numpy_sugar.linalg import cho_solve, ddot, dotd
+        from numpy_sugar.linalg import ddot, dotd, lu_solve
 
         self._update()
 
-        L = self._posterior.L()
+        LU = self._posterior.LU()
         Q = self._posterior.cov["QS"][0][0]
         ttau = self._site.tau
         teta = self._site.eta
 
         diff = teta - ttau * self._posterior.mean
 
         v0 = dot(dQS[0][0], dQS[1] * dot(dQS[0][0].T, diff))
-        v1 = ttau * dot(Q, cho_solve(L, dot(Q.T, diff)))
+        v1 = ttau * dot(Q, lu_solve(LU, dot(Q.T, diff)))
         dlml = 0.5 * dot(diff, v0)
         dlml -= dot(v0, v1)
         dlml += 0.5 * dot(v1, dot(dQS[0][0], dQS[1] * dot(dQS[0][0].T, v1)))
         dqs = ddot(dQS[1], dQS[0][0].T, left=True)
         diag = dotd(dQS[0][0], dqs)
         dlml -= 0.5 * sum(ttau * diag)
 
-        tmp = cho_solve(L, dot(ddot(Q.T, ttau, left=False), dQS[0][0]))
+        tmp = lu_solve(LU, dot(ddot(Q.T, ttau, left=False), dQS[0][0]))
         dlml += 0.5 * sum(ttau * dotd(Q, dot(tmp, dqs)))
 
         return dlml
 
     def lml_derivative_over_mean(self, dm):
-        from numpy_sugar.linalg import cho_solve
+        from numpy_sugar.linalg import lu_solve
 
         self._update()
 
-        L = self._posterior.L()
+        LU = self._posterior.LU()
         Q = self._posterior.cov["QS"][0][0]
         ttau = self._site.tau
         teta = self._site.eta
 
         diff = teta - ttau * self._posterior.mean
 
         dlml = dot(diff, dm)
-        dlml -= dot(diff, dot(Q, cho_solve(L, dot(Q.T, (ttau * dm.T).T))))
+        dlml -= dot(diff, dot(Q, lu_solve(LU, dot(Q.T, (ttau * dm.T).T))))
 
         return dlml
 

glimix_core/_ep/linear_kernel.py

@@ -25,12 +25,14 @@ def __init__(self, nsites, rtol=1.49e-05, atol=1.49e-08):
         )
 
     def lml(self):
+        from numpy_sugar.linalg import lu_slogdet
+
         if self._cache["lml"] is not None:
             return self._cache["lml"]
 
         self._update()
 
-        L = self._posterior.L()
+        LU = self._posterior.LU()
         LQt = self._posterior.LQt()
         cov = self._posterior.cov
         Q = cov["QS"][0][0]
@@ -49,7 +51,7 @@ def lml(self):
         dif = 2 * A * teta - A * ttau * m
 
         lml = [
-            -log(L.diagonal()).sum(),
+            -lu_slogdet(LU)[1],
             -0.5 * sum(log(s * S)),
             +0.5 * sum(log(A)),
             # lml += 0.5 * sum(log(ttau)),

glimix_core/_ep/posterior.py

@@ -1,4 +1,5 @@
-from numpy import dot, empty, empty_like, sqrt, sum as npsum
+from numpy import dot, empty, empty_like, sqrt
+from numpy import sum as npsum
 
 
 class Posterior(object):
@@ -23,7 +24,7 @@ def __init__(self, site):
         self._NxR_data = None
         self._RxN_data = None
         self._RxR_data = None
-        self._L_cache = None
+        self._LU_cache = None
         self._LQt_cache = None
         self._QS_cache = None
         self._AQ_cache = None
@@ -53,7 +54,7 @@ def _RxR(self):
         return self._RxR_data
 
     def _flush_cache(self):
-        self._L_cache = None
+        self._LU_cache = None
         self._LQt_cache = None
         self._QS_cache = None
         self._AQ_cache = None
@@ -104,37 +105,37 @@ def cov(self, v):
         self._initialize()
         self._cov = v
 
-    def L(self):
-        r"""Cholesky decomposition of :math:`\mathrm B`.
+    def LU(self):
+        r"""LU factor of :math:`\mathrm B`.
 
         .. math::
 
             \mathrm B = \mathrm Q^{\intercal}\tilde{\mathrm{T}}\mathrm Q
                 + \mathrm{S}^{-1}
         """
         from numpy_sugar.linalg import ddot, sum2diag
-        from scipy.linalg import cho_factor
+        from scipy.linalg import lu_factor
 
-        if self._L_cache is not None:
-            return self._L_cache
+        if self._LU_cache is not None:
+            return self._LU_cache
 
         Q = self._cov["QS"][0][0]
         S = self._cov["QS"][1]
         B = dot(Q.T, ddot(self._site.tau, Q, left=True))
         sum2diag(B, 1.0 / S, out=B)
-        self._L_cache = cho_factor(B, lower=True)[0]
-        return self._L_cache
+        self._LU_cache = lu_factor(B, overwrite_a=True, check_finite=False)
+        return self._LU_cache
 
     def LQt(self):
-        from numpy_sugar.linalg import cho_solve
+        from numpy_sugar.linalg import lu_solve
 
         if self._LQt_cache is not None:
             return self._LQt_cache
 
         L = self.L()
         Q = self._cov["QS"][0][0]
 
-        self._LQt_cache = cho_solve(L, Q.T)
+        self._LQt_cache = lu_solve(L, Q.T)
         return self._LQt_cache
 
     def QS(self):

glimix_core/_ep/posterior_linear_kernel.py

@@ -3,13 +3,6 @@
 from .posterior import Posterior
 
 
-def _cho_factor(B):
-    from scipy.linalg import cho_factor
-
-    B = cho_factor(B, overwrite_a=True, lower=True, check_finite=False)[0]
-    return B
-
-
 class PosteriorLinearKernel(Posterior):
     """
     EP posterior.
@@ -62,18 +55,19 @@ def QSQtATQLQtA(self):
         self._QSQtATQLQtA_cache = dotd(QS, dot(dot(ddot(AQ.T, T), Q), LQtA))
         return self._QSQtATQLQtA_cache
 
-    def L(self):
-        r"""Cholesky decomposition of :math:`\mathrm B`.
+    def LU(self):
+        r"""LU factor of :math:`\mathrm B`.
 
         .. math::
 
             \mathrm B = \mathrm Q^{\intercal}\tilde{\mathrm{T}}\mathrm Q
                 + \mathrm{S}^{-1}
         """
         from numpy_sugar.linalg import ddot, sum2diag
+        from scipy.linalg import lu_factor
 
-        if self._L_cache is not None:
-            return self._L_cache
+        if self._LU_cache is not None:
+            return self._LU_cache
 
         s = self._cov["scale"]
         d = self._cov["delta"]
@@ -84,8 +78,8 @@ def L(self):
         B = dot(Q.T, self._NxR, out=self._RxR)
         B *= 1 - d
         sum2diag(B, 1.0 / S / s, out=B)
-        self._L_cache = _cho_factor(B)
-        return self._L_cache
+        self._LU_cache = lu_factor(B, overwrite_a=True, check_finite=False)
+        return self._LU_cache
 
     def update(self):
         from numpy_sugar.linalg import ddot, dotd

glimix_core/glmm/test/test_glmm_expfam.py

@@ -16,7 +16,6 @@
 from numpy.random import RandomState
 from numpy.testing import assert_, assert_allclose
 from numpy_sugar.linalg import economic_qs, economic_qs_linear
-from pandas import DataFrame
 
 from glimix_core.example import linear_eye_cov
 from glimix_core.glmm import GLMMExpFam, GLMMNormal
@@ -606,9 +605,6 @@ def test_glmmexpfam_poisson():
     offset = ones(n) * random.randn()
     age = random.randint(16, 75, n)
     M = stack((offset, age), axis=1)
-    M = DataFrame(stack([offset, age], axis=1), columns=["offset", "age"])
-    M["sample"] = [f"sample{i}" for i in range(n)]
-    M = M.set_index("sample")
 
     # genetic variants
     G = random.randn(n, 4)

glimix_core/lmm/_kron2sum.py

@@ -196,7 +196,7 @@ def C0(self):
         C0 : ndarray
             C₀.
         """
-        return self._cov.C0.value() / (self._G_norm ** 2)
+        return self._cov.C0.value() / (self._G_norm**2)
 
     @property
     def C1(self):
@@ -419,8 +419,8 @@ def _XY(self):
 
     @property
     def _terms(self):
-        from numpy_sugar.linalg import ddot, sum2diag
-        from scipy.linalg import cho_factor, cho_solve
+        from numpy_sugar.linalg import ddot, lu_solve, sum2diag
+        from scipy.linalg import lu_factor
 
         if self._cache["terms"] is not None:
             return self._cache["terms"]
@@ -439,7 +439,7 @@ def _terms(self):
 
         Z = kron(L0.T @ WL0, self._GG)
         Z = sum2diag(Z, 1)
-        Lz = cho_factor(Z, lower=True)
+        Lz = lu_factor(Z, check_finite=False)
 
         # 𝐲ᵀR⁻¹𝐲 = vec(YW)ᵀ𝐲
         yRiy = (YW * self._Y).sum()
@@ -452,17 +452,17 @@ def _terms(self):
         # MᵀR⁻¹𝐲 = vec(XᵀYWA)
         MRiy = vec(self._XY @ WA)
 
-        ZiXRiM = cho_solve(Lz, XRiM)
-        ZiXRiy = cho_solve(Lz, XRiy)
+        ZiXRiM = lu_solve(Lz, XRiM)
+        ZiXRiy = lu_solve(Lz, XRiy)
 
         MRiXZiXRiy = ZiXRiM.T @ XRiy
         MRiXZiXRiM = XRiM.T @ ZiXRiM
 
         yKiy = yRiy - XRiy @ ZiXRiy
         MKiy = MRiy - MRiXZiXRiy
         H = MRiM - MRiXZiXRiM
-        Lh = cho_factor(H)
-        b = cho_solve(Lh, MKiy)
+        Lh = lu_factor(H, check_finite=False)
+        b = lu_solve(Lh, MKiy)
         B = unvec(b, (self.ncovariates, -1))
         self._mean.B = B
         XRim = XRiM @ b
@@ -571,7 +571,7 @@ def _lml_gradient(self):
         C1.Lu : ndarray
             Gradient of the log of the marginal likelihood over C₁ parameters.
         """
-        from scipy.linalg import cho_solve
+        from numpy_sugar.linalg import lu_solve
 
         terms = self._terms
         dC0 = self._cov.C0.gradient()["Lu"]
@@ -621,10 +621,10 @@ def _lml_gradient(self):
         yR0y = vec(_mdot(self._GY, WdC0)).T @ vec(self._GY @ W)
         yR1y = (YW.T * _mdot(self._Y, WdC1).T).T.sum(axis=(0, 1))
 
-        ZiXR0X = cho_solve(Lz, XR0X)
-        ZiXR1X = cho_solve(Lz, XR1X)
-        ZiXR0y = cho_solve(Lz, XR0y)
-        ZiXR1y = cho_solve(Lz, XR1y)
+        ZiXR0X = lu_solve(Lz, XR0X)
+        ZiXR1X = lu_solve(Lz, XR1X)
+        ZiXR0y = lu_solve(Lz, XR0y)
+        ZiXR1y = lu_solve(Lz, XR1y)
 
         # Mᵀ𝕂y = Mᵀ𝓡𝐲 - (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡𝐲) - (Mᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲)
         #       + (MᵀR⁻¹X)Z⁻¹(Xᵀ𝓡X)Z⁻¹(XᵀR⁻¹𝐲)
@@ -655,16 +655,16 @@ def _lml_gradient(self):
         XRim = XRiM @ b
         MRim = MRiM @ b
 
-        db = {"C0.Lu": cho_solve(Lh, MK0m - MK0y), "C1.Lu": cho_solve(Lh, MK1m - MK1y)}
+        db = {"C0.Lu": lu_solve(Lh, MK0m - MK0y), "C1.Lu": lu_solve(Lh, MK1m - MK1y)}
 
         grad = {
             "C0.Lu": -trace(WdC0) * self._trGG + trace(ZiXR0X),
             "C1.Lu": -trace(WdC1) * self.nsamples + trace(ZiXR1X),
         }
 
         if self._restricted:
-            grad["C0.Lu"] += cho_solve(Lh, MK0M).diagonal().sum(1)
-            grad["C1.Lu"] += cho_solve(Lh, MK1M).diagonal().sum(1)
+            grad["C0.Lu"] += lu_solve(Lh, MK0M).diagonal().sum(1)
+            grad["C1.Lu"] += lu_solve(Lh, MK1M).diagonal().sum(1)
 
         mKiM = MRim.T - XRim.T @ ZiXRiM
         yKiM = MRiy.T - XRiy.T @ ZiXRiM