mathmat.factor

  1"""MathMat Matrix Factorization Toolbox."""
  2
  3from .mathmat import MatrixTypeException, Matrix, DiagonalMatrix, Identity
  4import numpy as np
  5from scipy.linalg import lu as scipy_lu
  6from scipy.sparse.linalg import svds as sparse_svd, \
  7    splu as sparse_lu
  8from scipy.sparse import csc_array
  9
 10
 11def cholesky(M: Matrix):
 12    """Compute the Cholesky factorization $M = L L^H$."""
 13    if not M.is_hermitian():
 14        raise MatrixTypeException("Cholesky factorization", "Hermitian")
 15    if not M.is_posdef():
 16        raise MatrixTypeException(
 17            "Cholesky factorization", "positive-definite")
 18    L = np.linalg.cholesky(M.entries)
 19    L = Matrix(csc_array(L))
 20    L.set("dense", Matrix(L))
 21    L.set("triangular_L", True)
 22    M.set("cholesky", L)
 23    M.set("plu", (Identity(M.nc), L, L.transpose()))
 24    return L
 25
 26
 27def diagonalization(M: Matrix):
 28    """Compute the eigenvalue decomposition $M = P D P^{-1}$. Densifies."""
 29    if not M.is_diagonalizable():
 30        raise MatrixTypeException("diagonalize", "diagonalizable")
 31    if M.is_sparse():
 32        M = M.to_dense()
 33    P = Matrix(M.eigenvectors())
 34    D = DiagonalMatrix(M.eigenvalues())
 35
 36    M.set("diagonalization", (D, P))
 37    return D, P
 38
 39
 40def lu(M: Matrix, densify=False):
 41    """Compute the pivoted LU factorization $M = P L U$.
 42    If the matrix is sparse, returns an object that can be used to solve
 43    linear systems."""
 44    if not M.is_square():
 45        raise MatrixTypeException("LU factorization", "square")
 46    if M.is_sparse():
 47        if densify:
 48            M = M.to_dense()
 49        else:
 50            fct = sparse_lu(M.entries)
 51            M.set("lu_solver", fct.solve)
 52            return fct
 53    P, L, U = scipy_lu(M.entries)
 54    P = Matrix(P)
 55    L = Matrix(csc_array(L))
 56    L.set("dense", Matrix(L))
 57    L.set("triangular_L", True)
 58    U = Matrix(csc_array(U))
 59    U.set("dense", Matrix(U))
 60    U.set("triangular_R", True)
 61
 62    P.set("unitary", True)
 63    P.set("inverse", lambda: U.conjugate().transpose())
 64    M.set("plu", (P, L, U))
 65    return P, L, U
 66
 67
 68def qr(M: Matrix):
 69    """Compute the QR factorization $M = QR$."""
 70    Q, R = np.linalg.qr(M.entries, mode="reduced")
 71    Q = Matrix(Q)
 72    R = Matrix(R)
 73    R.set("triangular_U", True)
 74    # Store known properties of Q and R
 75    if M.is_complex():
 76        Q._computed.update({
 77            "unitary": True,
 78            "inverse": lambda: Q.conjugate().transpose()
 79        })
 80    else:
 81        Q._computed.update({
 82            "orthogonal": True,
 83            "inverse": lambda: Q.transpose()
 84        })
 85
 86    M.set("qr", (Q, R))
 87    return Q, R
 88
 89
 90def svd(M: Matrix, tolerance=1e-10, sp_k="max"):
 91    """Compute the SVD $M = U \\Sigma V^H$.
 92
 93    Note that for sparse matrices, there is a limit set by sparse.svds
 94    as to how many singular values can be found! `max` is set to half
 95    of the minimum dimension"""
 96    if M.is_sparse():
 97        # Use sparse SVD
 98        U, sigmas, V = sparse_svd(
 99            M.entries,
100            k=max(2, int(0.1*min(M.nr, M.nc)) if sp_k == "max" else sp_k),
101            solver="lobpcg")
102        order = np.flip(np.argsort(sigmas))
103        U = U[:, order]
104        sigmas = sigmas[order]
105        V = V[order, :]
106    else:
107        # Use dense SVD
108        U, sigmas, V = np.linalg.svd(M.entries)
109    U = Matrix(U)
110    if M.is_square():
111        # Use a sparse diagonal
112        S = DiagonalMatrix(sigmas)
113    else:
114        # Place the sigmas into the right shape
115        vals = np.zeros(shape=M.size)
116        vals[:len(sigmas), :len(sigmas)] = np.diag(sigmas)
117        S = Matrix(vals)
118    V = Matrix(V)
119    # Store the known properties of U, S, and V.
120    if M.is_complex():
121        U._computed.update({
122            "unitary": True,
123            "inverse": lambda: U.conjugate().transpose()
124        })
125        V._computed.update({
126            "unitary": True,
127            "inverse": lambda: V.conjugate().transpose()
128        })
129    else:
130        U._computed.update({
131            "orthogonal": True,
132            "inverse": lambda: U.transpose()
133        })
134        V._computed.update({
135            "orthogonal": True,
136            "inverse": lambda: V.transpose()
137        })
138
139    M.set("svd", (U, S, V))
140    M.set("sigmas", sigmas)
141    return U, S, V