First commit

.gitignore

@@ -0,0 +1,11 @@
+
+.DS_Store
+dist/egrssmatrices-0.0.1-py3-none-any.whl
+dist/egrssmatrices-0.0.1.tar.gz
+egrssmatrices.egg-info/SOURCES.txt
+egrssmatrices.egg-info/PKG-INFO
+egrssmatrices.egg-info/top_level.txt
+egrssmatrices/__pycache__
+egrssmatrices.egg-info/dependency_links.txt
+egrssmatrices/ideas.py
+tests/__pycache__

README.md

@@ -1,2 +1,19 @@
 # egrssmatrices
-
+
+## Description
+A package for efficiently computing with symmetric extended generator representable semiseparable matrices (EGRSS Matrices) and a variant with added diagonal terms. In short this means matrices of the form
+
+$$
+K = \text{\textbf{tril}}(UV^\top) + \text{\textbf{triu}}\left(VU^\top,1\right), \quad U,V\in\mathbb{R}^{p\times n}
+$$
+
+$$
+K = \text{\textbf{tril}}(UV^\top) + \text{\textbf{triu}}\left(VU^\top,1\right) + \text{\textbf{diag}}(d), \quad U,V\in\mathbb{R}^{p\times n},\ d\in\mathbb{R}^n
+$$
+
+All implemented algorithms (multiplication, Cholesky factorization, forward/backward substitution as well as various traces and determinants) scales with $O(p^kn)$. Since $p \ll n$ this result in very scalable computations.
+
+A more in-depth descriptions of the algorithms can be found in [1] or [here](https://github.com/mipals/SmoothingSplines.jl/blob/master/mt_mikkel_paltorp.pdf).
+
+## References
+[1] M. S. Andersen and T. Chen, “Smoothing Splines and Rank Structured Matrices: Revisiting the Spline Kernel,” SIAM Journal on Matrix Analysis and Applications, 2020.

egrssmatrices/egrssmatrix.py

@@ -0,0 +1,158 @@
+import numpy as np
+
+class egrssmatrix:
+    def __init__(self,Ut,Vt,d=None):
+        self.Ut = Ut
+        self.Vt = Vt
+        self.p, self.n = Ut.shape
+        self.d  = d
+        self.Wt = None
+        self.c  = None
+        self.Yt = None
+        self.Zt = None
+        self.shape = (self.n, self.n)
+        self.dtype = Ut.dtype
+
+    def __repr__(self):
+        return f"egrssmatrix(p={self.p},n={self.n})"
+
+    def __add__(self,other):
+        return egrssmatrix(np.hstack((self.Ut,other.Ut)), np.hstack((self.Ut,other.Ut)))
+
+    def __matmul__(self,other):
+        assert self.n == other.shape[0], f"dimension error. Expected length to be {self.n} got {other.shape[0]}"
+        ubar = self.Ut @ other
+        vbar = np.zeros(ubar.shape)
+        yk = np.zeros(other.shape)
+        for k in range(self.n):
+            vbar = vbar + self.Vt[:,k]*other[k]
+            ubar = ubar - self.Ut[:,k]*other[k]
+            yk[k] = np.dot(self.Ut[:,k],vbar) + np.dot(self.Vt[:,k],ubar)
+        if self.d is not None:
+            yk = yk + self.d*other
+        return yk
+
+    def __getitem__(self, key):
+        if key[0] > key[1]:
+            return np.dot(self.Ut[:,key[0]], self.Vt[:,key[1]])
+        elif key[0] == key[1]:
+            if self.d is not None:
+                return np.dot(self.Ut[:,key[0]], self.Vt[:,key[1]]) + self.d[key[0]]
+            else:
+                return np.dot(self.Ut[:,key[0]], self.Vt[:,key[1]])
+        return np.dot(self.Vt[:,key[0]], self.Ut[:,key[1]])
+
+    def __setitem__(self,key,value):
+        raise IndexError('You can not set index value of an egrssmatrix')
+
+    def matvec(self,other):
+        return self@other
+
+    def full(self):
+        K = np.tril(self.Ut.T@self.Vt,-1)
+        K = K + K.T + np.diag(np.sum(self.Ut*self.Vt,axis=0))
+        if self.d is not None:
+            if np.isscalar(self.d):
+                K = K + np.diag(self.d*np.ones(self.n))
+            else: 
+                K = K + np.diag(self.d)
+        return K
+
+    def cholesky(self):
+        wbar = np.zeros((self.p,self.p))
+        Wt = self.Vt.copy()
+        if self.d is None:
+            for k in range(self.n):
+                Wt[:,k] = Wt[:,k] - wbar@self.Ut[:,k]
+                Wt[:,k] = Wt[:,k]/np.sqrt(np.dot(self.Ut[:,k],Wt[:,k]))
+                wbar   = wbar + np.outer(Wt[:,k],Wt[:,k])
+            self.Wt = Wt
+            self.c = np.sum(self.Ut*self.Wt,axis=0)
+        else:
+            if np.isscalar(self.d):
+                c = self.d*np.ones(self.n)
+            else:
+                c = self.d.copy()
+
+            for k in range(self.n):
+                Wt[:,k] = Wt[:,k] - wbar@self.Ut[:,k]
+                c[k] = np.sqrt(np.dot(self.Ut[:,k],Wt[:,k]) + c[k])
+                Wt[:,k] = Wt[:,k]/c[k]
+                wbar   = wbar + np.outer(Wt[:,k],Wt[:,k])
+            self.Wt = Wt
+            self.c = c
+
+    def solve(self,other):
+        return self.backward(self.forward(other))
+
+    def forward(self,other):
+        if self.Wt is None:
+            self.cholesky()
+        wbar = np.zeros(self.p)
+        x = np.zeros(self.n)
+        for k in range(self.n):
+            x[k] = (other[k] - np.dot(self.Ut[:,k],wbar))/self.c[k]
+            wbar = wbar + self.Wt[:,k]*x[k]
+        return x
+
+    def backward(self,other):
+        if self.Wt is None:
+            self.cholesky()
+        ubar = np.zeros(self.p)
+        x = np.zeros(self.n)
+        for k in reversed(range(self.n)):
+            x[k] = (other[k] - np.dot(self.Wt[:,k],ubar))/self.c[k]
+            ubar = ubar + self.Ut[:,k]*x[k]
+        return x
+
+    def cholesky_inverse(self):
+        if self.Yt is not None and self.Zt is not None:
+            return
+        if self.Wt is None:
+            self.cholesky()
+        Yt = self.Ut.copy()
+        Zt = self.Wt.copy()
+
+        for k in range(self.p):
+            Yt[k,:] = self.forward(self.Ut[k,:])
+            Zt[k,:] = self.backward(self.Wt[k,:])
+
+        Zt = np.linalg.solve(Zt@self.Ut.T - np.eye(self.p), Zt)
+
+        self.Yt = Yt
+        self.Zt = Zt
+
+    def logdet(self):
+        if self.c is None:
+            self.cholesky()
+        return 2*np.sum(np.log(self.c))
+
+    def trace_inverse(self):
+        if self.Yt is None and self.Zt is None:
+            self.cholesky_inverse()
+
+        Ybar = np.zeros((self.p,self.p))
+        b = 0
+        for k in reversed(range(self.n)):
+            b = b + self.c[k]**(-2) + self.Zt[:,k]@(Ybar@self.Zt[:,k])
+            Ybar = Ybar + np.outer(self.Yt[:,k],self.Yt[:,k])
+
+        return b
+
+    def trace_inverse_product(self):
+        if self.Yt is None and self.Zt is None:
+            self.cholesky_inverse()
+        b = 0
+        P = np.zeros((self.p,self.p))
+        R = np.zeros((self.p,self.p))
+
+        for k in range(self.n):
+            b = b + np.dot(self.Yt[:,k],P@self.Yt[:,k]) + \
+                  2*np.dot(self.Yt[:,k],R@self.Ut[:,k])/self.c[k] + \
+                    np.dot(self.Ut[:,k],self.Vt[:,k])/(self.c[k]**2)
+            P = P + np.dot(self.Ut[:,k],self.Vt[:,k])*np.outer(self.Zt[:,k],self.Zt[:,k]) + \
+                    np.outer(self.Zt[:,k],R@self.Ut[:,k]) + \
+                    np.outer(R@self.Ut[:,k],self.Zt[:,k])
+            R = R + np.outer(self.Zt[:,k],self.Vt[:,k])
+
+        return b

egrssmatrices/spline_kernel.py

@@ -0,0 +1,7 @@
+import numpy as np
+
+def spline_kernel(t,p):
+    T = np.vander(t, 2*p)/np.flip(np.cumprod(np.hstack([1,np.arange(1,2*p)])))
+    Ut = T[:,p:].T
+    Vt = (np.fliplr(T[:,:p])*((-1)**np.arange(0,p))).T
+    return Ut,Vt

pyproject.toml

@@ -0,0 +1,21 @@
+[build-system]
+requires = ["hatchling", "numpy"]
+build-backend = "hatchling.build"
+
+[project]
+name = "egrssmatrices"
+version = "0.0.1"
+authors = [
+  { name="Mikkel Paltorp", email="mikkel.paltorp@gmail.com" },
+]
+description = "A small example package"
+readme = "README.md"
+requires-python = ">=3.8"
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "License :: OSI Approved :: MIT License",
+    "Operating System :: OS Independent",
+]
+
+[project.urls]
+"Homepage" = "https://github.com/mipals/egrssmatrices"

tests/test_egrssmatrix.py

@@ -0,0 +1,107 @@
+import unittest
+import numpy as np
+import numpy.testing as npt
+from egrssmatrices.egrssmatrix import egrssmatrix
+from egrssmatrices.spline_kernel import spline_kernel
+
+class testegrssmatrices(unittest.TestCase):
+
+    def test_matmul(self):
+        # Setting up the system
+        p,n = 2,50
+        t = np.linspace(0.02,1.0,n)
+        Ut,Vt = spline_kernel(t,p)
+        M  = egrssmatrix(Ut,Vt)
+        Md = egrssmatrix(Ut,Vt,t) # Vector diagonal   
+        Ms = egrssmatrix(Ut,Vt,1) # Scalar diagonal
+        Mfull  = M.full()
+        Mdfull = Md.full()
+        Msfull = Ms.full()
+        # Testing Multiplcation
+        x = np.ones(n)
+        y  = M@x
+        yd = Md@x
+        ys = Ms@x
+        yfull  = Mfull@x
+        ydfull = Mdfull@x
+        ysfull = Msfull@x
+        npt.assert_almost_equal(y,yfull)
+        npt.assert_almost_equal(yd,ydfull)
+        npt.assert_almost_equal(ys,ysfull)
+        # Testing __getitem__
+        npt.assert_equal(Mfull[1,10],  M[1,10])
+        npt.assert_equal(Mfull[10,10], M[10,10])
+        npt.assert_equal(Mfull[10,1],  M[10,1])
+        npt.assert_equal(Mdfull[1,10], Md[1,10])
+        npt.assert_equal(Mdfull[10,10],Md[10,10])
+        npt.assert_equal(Mdfull[10,1], Md[10,1])
+
+    def test_cholesky(self):
+        # Setting up the system
+        p,n = 2,50
+        t = np.linspace(0.02,1.0,n)
+        Ut,Vt = spline_kernel(t,p)
+        M  = egrssmatrix(Ut,Vt)
+        Md = egrssmatrix(Ut,Vt,t)
+        Ms = egrssmatrix(Ut,Vt,1)
+        Mfull  = M.full()
+        Mdfull = Md.full()
+        Msfull = Ms.full()
+        # Setting up right-hand sides
+        x  = np.ones(n)
+        y  = M@x
+        yd = Md@x
+        ys = Ms@x
+        # Solving the linear systems
+        b  = M.solve(y)
+        bd = Md.solve(yd)
+        bs = Ms.solve(ys)
+        bfull  = np.linalg.solve(Mfull,y)
+        bdfull = np.linalg.solve(Mdfull,yd)
+        bsfull = np.linalg.solve(Msfull,ys)
+        npt.assert_almost_equal(b,bfull)
+        npt.assert_almost_equal(bd,bdfull)
+        npt.assert_almost_equal(bs,bsfull)
+        # Testing traces
+        dtracefull = np.matrix.trace(np.linalg.solve(Mdfull,Mfull))
+        stracefull = np.matrix.trace(np.linalg.solve(Msfull,Mfull))
+        dtrace = Md.trace_inverse_product()
+        strace = Ms.trace_inverse_product()
+        npt.assert_almost_equal(dtrace,dtracefull)
+        npt.assert_almost_equal(strace,stracefull)
+
+    def test_cholesky_inverse(self):
+        # Setting up the system
+        p,n = 2,50
+        t = np.linspace(0.02,1.0,n)
+        Ut,Vt = spline_kernel(t,p)
+        Md = egrssmatrix(Ut,Vt,t)
+        Ms = egrssmatrix(Ut,Vt,1)
+        Mdfull = Md.full()
+        Msfull = Ms.full()
+        # Computing the inverse of the Cholesky factor
+        Md.cholesky_inverse()
+        Ms.cholesky_inverse()
+        Tdinv = np.tril(Md.Yt.T@Md.Zt,-1) + np.diag(1/Md.c)
+        Tsinv = np.tril(Ms.Yt.T@Ms.Zt,-1) + np.diag(1/Ms.c)
+        npt.assert_almost_equal(Tdinv,np.linalg.inv(np.linalg.cholesky(Mdfull)))
+        npt.assert_almost_equal(Tsinv,np.linalg.inv(np.linalg.cholesky(Msfull)))
+
+    def test_logdet(self):
+        # Setting up the system
+        p,n = 2,50
+        t = np.linspace(0.02,1.0,n)
+        Ut,Vt = spline_kernel(t,p)
+        M  = egrssmatrix(Ut,Vt)
+        Md = egrssmatrix(Ut,Vt,t)
+        Ms = egrssmatrix(Ut,Vt,1)
+        Mfull  = M.full()
+        Mdfull = Md.full()
+        Msfull = Ms.full()
+        # Computing the log-determinant
+        npt.assert_almost_equal(M.logdet(),np.linalg.slogdet(Mfull)[1])
+        npt.assert_almost_equal(Md.logdet(),np.linalg.slogdet(Mdfull)[1])
+        npt.assert_almost_equal(Ms.logdet(),np.linalg.slogdet(Msfull)[1])
+
+if __name__ == '__main__':
+    unittest.main()