mathmat.mathmat

  1"""MathMat is a mathematically aware matrix library for Python.
  2
  3See the Matrix class for more details.
  4"""
  5
  6import lazy_import
  7from functools import wraps
  8from sys import float_info
  9from scipy.sparse.linalg import eigs as sparse_eigs, eigsh as sparse_eigsh, \
 10    norm as sparse_norm
 11from scipy.sparse import identity, issparse, dia_array, csc_array, \
 12    tril as sparse_tril, triu as sparse_triu
 13from scipy.linalg import issymmetric, ishermitian, eig, eigh, eigvals, \
 14    eigvalsh, inv
 15import numpy as np
 16
 17
 18factor = lazy_import.lazy_module("mathmat.factor")
 19solve = lazy_import.lazy_module("mathmat.solve")
 20
 21EPS = float_info.epsilon
 22
 23
 24def uniquetol(arr, tolerance):
 25    """Find and return the unique elements of `arr` within a `tolerance`."""
 26    unique = np.unique(arr)
 27    if len(unique) <= 1:
 28        return unique
 29    u_diff = np.abs(np.diff(unique))
 30    where = np.argwhere(u_diff > tolerance)
 31    return np.insert(unique[where + 1], 0, unique[0])
 32
 33
 34class MatrixSizeException(Exception):
 35    """A `MatrixSizeException` is raised when two matrices have
 36    incompatible sizes for an operation, usually multiplication."""
 37
 38    def __init__(self, A, B, operation):
 39        super().__init__(
 40            "Incompatible matrix sizes {}x{} and {}x{} for {}.".format(
 41                A.nr, A.nc, B.nr, B.nc, operation))
 42
 43
 44class MatrixTypeException(Exception):
 45    """A `MatrixTypeException` is raised when a matrix does not satisfy
 46    the requirements of an attempted operation."""
 47
 48    def __init__(self, operation, req_type):
 49        super().__init__(
 50            "The matrix must be {} to perform: {}.".format(
 51                req_type, operation))
 52
 53
 54def _computable_property(name):
 55    """Decorate a function as defining a Computable Property of a Matrix."""
 56    def decorator(func):
 57        @wraps(func)
 58        def property_wrapper(*args, **kwargs):
 59            """Check if the Property already has a value."""
 60            if name not in args[0]._computed:
 61                # No value exists, perform the computation.
 62                args[0].set(name, func(*args, **kwargs))
 63            value = args[0]._computed.get(name, None)
 64            if callable(value) and not name == "lin_solver":
 65                # A value can be generated using another function.
 66                value = value()
 67                args[0].set(name, value)
 68            return value
 69        return property_wrapper
 70    return decorator
 71
 72
 73class Matrix:
 74    """A Matrix object is an immutable representation of a matrix.
 75
 76    In addition to the entries, the Matrix keeps track of
 77    computations which have been performed on it, and reuses
 78    existing computations whenever possible. Thus, a Matrix
 79    object is "aware" of linear algebra theorems.
 80
 81    Each Computable Property can be either:
 82     1. unset (default), in which case the corresponding
 83        function will be called to compute it.
 84     2. set to a callable, in which case that callable
 85        will be called to compute it.
 86     3. set to a value, which will be returned immediately.
 87
 88    List of Computable Properties:
 89     - complex        : `True` if $M$ has complex entries
 90     - conjugate      : The complex conjugate of $M$
 91     - cond_2         : Condition number in 2-norm
 92     - determinant    : `det(M)`
 93     - diagonal       : `True` if $M$ is zero except for main diag.
 94     - diagonalizable : `True` if $M = P D P^{-1}$
 95     - eigenvalues    : The eigenvalues of $M$
 96     - eigenvectors   : The right eigenvectors of M
 97     - hermitian      : `True` if conjugate transpose of $M$ equals $M$
 98     - inverse        : if it exists, M^{-1}
 99     - invertible     : `True` if M^{-1} exists
100     - lin_solver     : A function to solve Ax = b for x.
101     - norm_2         : The 2-norm of the matrix
102     - normal         : `True` if $M M^H = M^H M$
103     - nullity        : `dim(Kernel(M))`
104     - orthogonal     : `True` if $M^T = M^{-1}$
105     - posdef         : `True` if $M$ is positive definite
106     - qr             : The `(Q, R)` factorization of $M$
107     - rank           : `dim(Col(M))`
108     - sigmas         : Singular values of $M$
109     - sparsity       : The fraction of zero entries in $M$
110     - square         : `True` if $M$ is square
111     - svd            : The SVD `(U, S, V)` of $M$
112     - symmetric      : `True` if $M = M^T$
113     - trace          : The trace of $M$
114     - transpose      : $M^T$
115     - triangular_L   : `True` if $M$ is lower triangular.
116     - triangular_U   : `True` if $M$ is upper triangular.
117     - unitary        : `True` if $M^H = M^{-1}$
118
119    Non-Mathematical:
120     - sparse         : `True` if the matrix is stored as a
121                        SciPy Sparse object.
122     - to_dense       : A dense representation of $M$
123     - to_sparse      : A sparse representation of $M$
124    """
125
126    def set(self, prop, value):
127        self._computed[prop] = value
128
129    def __init__(self, arr):
130        """Initialize a Matrix object.
131
132        The argument `arr` can be a Numpy array, Scipy sparse matrix, list of
133        Vector objects, or an object convertible to a Numpy array.
134
135        If the resulting array is of dimension 1 in either axis, the
136        Matrix automatically becomes a Vector.
137        """
138        if isinstance(arr, Matrix):
139            arr = arr.entries
140        if issparse(arr):
141            # This is a SciPy sparse object. Preserve it.
142            self.entries = arr
143        elif (type(arr) is list or type(arr) is tuple) and len(arr) > 0 and \
144                isinstance(arr[0], Vector):
145            # This is a list of Vectors. Join them into a Matrix.
146            any_complex = sum(v.is_complex() for v in arr) > 0
147            if arr[0].nr == 1:
148                # Treat all as row vectors
149                self.entries = np.empty(
150                    shape=(len(arr), len(arr[0])),
151                    dtype=np.complex_ if any_complex else np.double)
152                for j in range(len(arr)):
153                    if len(arr[j]) != len(arr[0]):
154                        raise MatrixSizeException(
155                            arr[j], arr[0], "initialize Matrix from rows")
156                    self.entries[j, :] = arr[j].entries.ravel()
157            else:
158                # Treat all as column vectors
159                self.entries = np.empty(
160                    shape=(len(arr[0]), len(arr)),
161                    dtype=np.complex_ if any_complex else np.double)
162                for j in range(len(arr)):
163                    if len(arr[j]) != len(arr[0]):
164                        raise MatrixSizeException(
165                            arr[j], arr[0], "initialize Matrix from rows")
166                    self.entries[:, j] = arr[j].entries.ravel()
167        else:
168            # This is another kind of object, hopefully numpy can handle it.
169            self.entries = np.array(arr)
170        if len(self.entries.shape) == 1:
171            # Reshape ambiguous arrays into column vectors.
172            self.entries = self.entries.reshape((len(self.entries), 1))
173        if min(self.entries.shape) == 1 and not isinstance(self, Vector):
174            # Automatically become a Vector if has a dimension of 1
175            self.__dict__ = Vector(arr).__dict__
176            self.__class__ = Vector
177
178        # Initialize a blank dictionary of Computed Properties
179        self._computed = {}
180
181    """Directly accessible properties of the matrix"""
182
183    @property
184    def size(self) -> tuple:
185        """The shape of the underlying matrix representation."""
186        return self.entries.shape
187
188    @property
189    def nr(self) -> int:
190        """The number of rows in the Matrix."""
191        return self.size[0]
192
193    @property
194    def nc(self) -> int:
195        """The number of columns in the Matrix."""
196        return self.size[1]
197
198    def __str__(self):
199        return str(self.entries)
200
201    def __repr__(self):
202        return "{}x{} Matrix ({} properties)".format(
203            self.nr, self.nc, len(self._computed))
204
205    """Computable Properties corresponding to "checks".
206    All correspond to matrix "types" and thus are queried
207    using a method called `is_<type>()`."""
208
209    @_computable_property("sparse")
210    def is_sparse(self):
211        """`True` if the underlying matrix representation is SciPy Sparse."""
212        return issparse(self.entries)
213
214    @_computable_property("complex")
215    def is_complex(self):
216        """`True` if the matrix has any non-zero complex entry."""
217        return np.iscomplex(self.entries).any()
218
219    @_computable_property("diagonal")
220    def is_diagonal(self, tolerance=1e-10):
221        """`True` if the non-diagonal elements of the matrix are within
222        `tolerance` of zero."""
223        check = False
224        if self._computed.get("triangular_L", False) and \
225                self._computed.get("triangular_U", False):
226            return True
227        if self.is_square():
228            if self.is_sparse():
229                check = (self.entries -
230                         dia_array(
231                             (self.entries.diagonal().T, 0),
232                             shape=(self.nr, self.nc))).count_nonzero() == 0
233            else:
234                # Fast check for square matrices, @Daniel F on Stackoverflow
235                test = self.entries.reshape(-1)[:-1].reshape(
236                    self.nr-1, self.nc+1)
237                check = ~np.any(test[:, 1:])
238        else:
239            check = np.allclose(
240                self.entries - np.diag(self.entries, 0), 0, atol=tolerance)
241        if check:
242            # Automatically convert to a DiagonalMatrix
243            self.entries = dia_array((self.diagonal(0).entries.T, 0),
244                                     shape=(self.nr, self.nc))
245            self.__class__ = DiagonalMatrix
246        return check
247
248    @_computable_property("diagonalizable")
249    def is_diagonalizable(self, tolerance=1e-10):
250        """`True` if the matrix can be diagonalized, $M = P D P^{-1}$.
251        Does not necessarily compute the diagonalization. Densifies."""
252        if not self.is_square():
253            return False
254        if self.is_sparse():
255            return self.to_dense().is_diagonalizable()
256        if "eigenvalues" in self._computed:
257            if len(uniquetol(self.eigenvalues(), tolerance)) == self.nc:
258                # N distinct eigenvalues is sufficient for diagonalizability.
259                return True
260        if self.is_normal():
261            # Normal matrices are diagonalizable.
262            return True
263        # Check the linear independence of the eigenvectors through the SVD
264        return Matrix(self.eigenvectors()).rank(tolerance) == self.nc
265
266    @_computable_property("hermitian")
267    def is_hermitian(self, tolerance=1e-10):
268        """`True` if the matrix is Hermitian, $M^H = M$."""
269        if not self.is_square():
270            return False
271        if self.is_complex():
272            # Call the SciPy ishermitian function
273            return ishermitian(self.entries, atol=tolerance)
274        # For real matrices, equivalent to symmetric
275        return self.is_symmetric()
276
277    @_computable_property("invertible")
278    def is_invertible(self, tolerance=1e-10):
279        """`True` if the matrix has an inverse $M^{-1} M = Id$.
280        Does not compute the inverse."""
281        if "inverse" in self._computed:
282            return True
283        if "sigmas" in self._computed:
284            # No zero singular values => invertible
285            return np.isclose(self.sigmas(), 0, atol=tolerance).sum() == 0
286        if (self.is_triangular_L() or self.is_triangular_U()) or \
287                (not self.is_sparse() and "eigenvalues" in self._computed):
288            # No zero eigenvalues => invertible
289            return np.isclose(self.eigenvalues(), 0, atol=tolerance).sum() == 0
290        # Square matrices with a finite condition number are invertible
291        return self.is_square() and self.cond_2() < 1 / EPS
292
293    @_computable_property("normal")
294    def is_normal(self, tolerance=1e-10):
295        """`True` if the matrix satisfies the normal equation $M^H M = M M^H$.
296        Computes the conjugate transpose."""
297        if not self.is_square():
298            return False
299        CT = self.conjugate().transpose()
300        return (self @ CT).equals((CT @ self), tolerance)
301
302    @_computable_property("orthogonal")
303    def is_orthogonal(self, tolerance=1e-10):
304        """`True` if the matrix is orthogonal, $M^T = M^{-1}$.
305        Computes the transpose and identifies the inverse if `True`."""
306        if not self.is_square():
307            raise MatrixTypeException("orthogonality check", "square")
308        prod = self @ self.transpose()
309        check = prod.equals(Identity(prod.nc), tolerance)
310
311        if check:
312            # The inverse is now known.
313            self.set("inverse", self.transpose())
314            if not self.is_complex():
315                self.set("unitary", True)
316            self.transpose().set("orthogonal", True)
317            self.transpose().set("inverse", self)
318        return check
319
320    @_computable_property("posdef")
321    def is_posdef(self):
322        """`True` if the matrix is positive-definite.
323        Computes the eigenvalues. Densifies."""
324
325        if self.is_sparse():
326            return self.to_dense().is_posdef()
327        return self.is_symmetric() and np.all(self.eigenvalues() > 0)
328
329    def is_square(self):
330        """`True` if the matrix is square."""
331        return self.nr == self.nc
332
333    @_computable_property("symmetric")
334    def is_symmetric(self, tolerance=1e-10):
335        """`True` if the matrix is symmetric, $M^T = M$.
336        Computes the transpose only if the matrix is sparse."""
337        if not self.is_square():
338            return False
339        else:
340            if self.is_sparse():
341                return (abs(self.entries
342                            - self.transpose().entries) > tolerance).nnz == 0
343            # Call SciPy issymmetric
344            return issymmetric(self.entries, atol=tolerance)
345
346    @_computable_property("triangular_L")
347    def is_triangular_L(self, tolerance=1e-10):
348        """`True` if the matrix is lower (left) triangular."""
349        if self.is_sparse():
350            return (abs(self.entries
351                        - sparse_tril(self.entries)) > tolerance).sum() == 0
352        return np.allclose(self.entries, np.tril(self.entries), atol=tolerance)
353
354    @_computable_property("triangular_U")
355    def is_triangular_U(self, tolerance=1e-10):
356        """`True` if the matrix is upper (right) triangular."""
357        if self.is_sparse():
358            return (abs(self.entries
359                        - sparse_triu(self.entries)) > tolerance).sum() == 0
360        return np.allclose(self.entries, np.triu(self.entries), atol=tolerance)
361
362    @_computable_property("unitary")
363    def is_unitary(self, tolerance=1e-10):
364        """`True` if the matrix is unitary, $M^H = M^{-1}$.
365        Computes the conjugate transpose and identifies the inverse if True."""
366        if not self.is_square():
367            raise MatrixTypeException("unitary check", "square")
368        if not self.is_complex():
369            return self.is_orthogonal()
370        CT = self.conjugate().transpose()
371        prod = self @ CT
372        check = prod.equals(Identity(prod.nc), tolerance)
373
374        if check:
375            # The inverse is now known.
376            self.set("inverse", CT)
377            CT.set("unitary", True)
378            CT.set("inverse", self)
379        return check
380
381    """Computable Properties of the matrix that are not checks."""
382
383    @_computable_property("cond_2")
384    def cond_2(self):
385        """Compute the condition number $\\kappa_2 for inversion of the matrix.
386        Computes the singular values."""
387        sigmas = self.sigmas()
388        return sigmas[0] / sigmas[-1]
389
390    @_computable_property("determinant")
391    def determinant(self):
392        """Compute the determinant of the matrix. Densifies.
393        Computes the eigenvalues of triangular matrices."""
394        if not self.is_square():
395            raise MatrixTypeException("determinant", "square")
396        if self.is_sparse():
397            # For now, we have to go to a dense representation
398            return self.to_dense().determinant()
399        else:
400            if "eigenvalues" in self._computed or \
401                    (self._computed.get("triangular_L", False) or
402                     self._computed.get("triangular_U", False)):
403                # Use the eigenvalues if available
404                return np.prod(self.eigenvalues())
405            return np.linalg.det(self.entries)
406
407    @_computable_property("eigenvalues")
408    def eigenvalues(self, tolerance=1e-10):
409        """Compute the eigenvalues of the matrix.
410        Does not compute the eigenvectors, unless the matrix is sparse.
411        Checks whether $M$ is Hermitian, unless the matrix is triangular."""
412        if not self.is_square():
413            raise MatrixTypeException("eigenvalues", "square")
414        if self.is_triangular_L(tolerance) or self.is_triangular_U(tolerance):
415            # Triangular matrix eigenvalues are on the diagonal.
416            return np.sort(self.diagonal(0))
417        if self.is_sparse():
418            # Use sparse algorithms
419            self.eigenvectors(tolerance)
420            return self.eigenvalues()
421        if self.is_hermitian():
422            # Use the specialized algorithm for Hermitian matrices.
423            return eigvalsh(self.entries, check_finite=False)
424        vals = eigvals(self.entries, check_finite=False)
425        if np.iscomplex(vals).any():
426            return np.sort_complex(vals)
427        return np.sort(vals)
428
429    @_computable_property("eigenvectors")
430    def eigenvectors(self, tolerance=1e-10, sparse_k="max"):
431        """Compute the eigenvectors of the matrix.
432        Computes the eigenvalues. Checks whether $M$ is Hermitian.
433
434        Note that if the matrix has a sparse representation, then
435        only `sparse_k` many eigenvalues/vectors can be computed!
436        The `max` is 90% of the dimension of the matrix.
437        """
438        if not self.is_square():
439            raise MatrixTypeException("eigenvectors", "square")
440        if self.is_sparse():
441            sparse_k = min(self.nr - 1, int(0.9 * self.nr))
442            # Use the sparse algorithms
443            if self.is_hermitian():
444                # Use the specialized algorithm for Hermitian matrices
445                vals, vecs = sparse_eigsh(self.entries, k=sparse_k)
446            else:
447                vals, vecs = sparse_eigs(self.entries, k=sparse_k)
448            if not np.iscomplex(vals).any():
449                order = np.argsort(vals)
450                # Report the eigenvalues and vectors in increasing order
451                vals = vals[order]
452                vecs = vecs[:, order]
453        else:
454            if self.is_hermitian():
455                # Use the specialized algorithm for Hermitian matrices.
456                vals, vecs = eigh(self.entries, check_finite=False)
457            else:
458                vals, vecs = eig(self.entries, check_finite=False)
459                if not np.iscomplex(vals).any():
460                    order = np.argsort(vals)
461                    # Report the eigenvalues and vectors in increasing order
462                    vals = vals[order]
463                    vecs = vecs[:, order]
464        self.set("eigenvalues", vals)
465        return [Vector(vecs[:, j]) for j in range(len(vals))]
466
467    @_computable_property("norm_2")
468    def norm_2(self):
469        """Compute the 2-norm of the matrix $\\lVert M \\rVert_2."""
470        if self._computed.get("orthogonal", False) or \
471                self._computed.get("unitary", False):
472            # Unitary and Orthogonal matrices have norm 1.
473            return 1.0
474        if "sigmas" in self._computed:
475            # The norm is the largest singular value.
476            return self.sigmas()[0]
477        if self.is_sparse():
478            return sparse_norm(self.entries, ord=2)
479        return np.linalg.norm(self.entries, ord=2)
480
481    @_computable_property("nullity")
482    def nullity(self, tolerance=1e-10):
483        """Compute the nullity as the number of columns minus the rank."""
484        return self.nc - self.rank()
485
486    @_computable_property("rank")
487    def rank(self, tolerance=1e-10):
488        """Compute the rank of the matrix.
489        Computes the SVD if no inverse or eigenvalues are known."""
490        if "invertible" in self._computed or "inverse" in self._computed:
491            # Invertible => full rank
492            return self.nc
493        if not self.is_sparse() and "eigenvalues" in self._computed:
494            # Number of nonzero eigenvalues
495            return np.logical_not(
496                np.isclose(self.eigenvalues(), 0, atol=tolerance)).sum()
497        # Number of nonzero signular values
498        _, S, _ = factor.svd(self)
499        return np.logical_not(
500            np.isclose(S.diagonal(0).entries, 0, atol=tolerance)).sum()
501
502    @_computable_property("sigmas")
503    def sigmas(self):
504        """Compute the singular values of the matrix.
505        Computes the SVD if not positive-definite."""
506        if "posdef" in self._computed and self.is_posdef():
507            return np.flip(np.sqrt(self.eigenvalues()))
508        factor.svd(self)
509        return self.sigmas()
510
511    @_computable_property("sparsity")
512    def sparsity(self, tolerance=1e-10):
513        """Compute the sparsity of the matrix as the fraction of zero
514        entries out of all entries in the matrix."""
515        if self.is_sparse():
516            # Use the sparse methods if available
517            return len(self.entries.nnz) / (self.nr * self.nc)
518        return np.isclose(
519            self.entries, 0, atol=tolerance).sum() / (self.nr * self.nc)
520
521    @_computable_property("trace")
522    def trace(self):
523        """Compute the trace of the matrix."""
524        if self.is_sparse():
525            return self.entries.trace()
526        if "eigenvalues" in self._computed:
527            return np.sum(self.eigenvalues())
528        return np.trace(self.entries)
529
530    """Getter functions for elements of the matrix."""
531
532    def row(self, number):
533        """Return the 1-indexed row."""
534        if number < 0:
535            return Vector(self.entries[number, :])
536        return Vector(self.entries[number-1, :])
537
538    def column(self, number):
539        """Return the 1-indexed column."""
540        if number < 0:
541            return Vector(self.entries[:, number])
542        return Vector(self.entries[:, number-1])
543
544    def diagonal(self, number):
545        """Return the k-th diagonal."""
546        if self.is_sparse():
547            return Vector(self.entries.diagonal(number))
548        return Vector(np.diagonal(self.entries, number))
549
550    def entry(self, row, col):
551        """Return the 1-indexed entry $M_{ij}$."""
552        return self.entries[row if row < 0 else row - 1,
553                            col if col < 0 else col-1]
554
555    def equals(self, M, tolerance=1e-10):
556        """`True` if the matrix equals another elementwise,
557        within a tolerance."""
558        if type(M) is int or type(M) is float:
559            return np.allclose(self.entries, M, atol=tolerance)
560        if self.is_sparse() or M.is_sparse():
561            # Use subtraction if one or both are sparse
562            return np.allclose(self.entries - M.entries, 0, atol=tolerance)
563        return np.allclose(self.entries, M.entries, atol=tolerance)
564
565    def __eq__(self, M):
566        return self.equals(M)
567
568    """Transformations of the matrix."""
569
570    @_computable_property("dense_mat")
571    def to_dense(self):
572        """Switch from a sparse to a dense repsresentation."""
573        if self.is_sparse():
574            DM = Matrix(self.entries.todense())
575            DM._computed = self._computed
576            # Remove properties that are best recomputed
577            for k in ("cond_2", "dense", "eigenvalues", "eigenvectors", "rank",
578                      "invertible", "inverse", "sigmas", "sparse", "svd",
579                      "transpose"):
580                DM._computed.pop(k, None)
581            DM.set("sparse_mat", self)
582            return DM
583        return self
584
585    @_computable_property("sparse_mat")
586    def to_sparse(self, tolerance=1e-10):
587        """Switch from a dense to a sparse representation."""
588        if self.is_sparse():
589            return self
590        if self.is_diagonal():
591            return DiagonalMatrix(self.diagonal(0))
592        return Matrix(csc_array(self.entries))
593
594    @_computable_property("conjugate")
595    def conjugate(self):
596        """Compute the complex conjugate of the matrix."""
597        if self.is_sparse():
598            Conj = Matrix(self.entries.conjugate())
599        else:
600            Conj = Matrix(np.conjugate(self.entries))
601        Conj.set("conjugate", self)
602        if "triangular_L" in self._computed:
603            Conj.set("triangular_L", self.is_triangular_L())
604        if "triangular_U" in self._computed:
605            Conj.set("triangular_U", self.is_triangular_U())
606        if "diagonal" in self._computed:
607            Conj.set("diagonal", self.is_diagonal())
608        return Conj
609
610    @_computable_property("inverse")
611    def inverse(self):
612        """Compute the inverse of the matrix. Checks invertibility."""
613        if not self.is_invertible():
614            raise MatrixTypeException("inverse", "invertible")
615        if "qr" in self._computed:
616            # Use the QR if it exists
617            Q, R = self.qr()
618            Inv = R.inverse() @ Q.transpose()
619        elif "svd" in self._computed:
620            # Use the SVD if it exists
621            U, S, V = self.svd()
622            Inv = V.conjugate().transpose() @ S.inverse() \
623                @ U.conjugate().transpose()
624        else:
625            # Nothing known, compute inverse
626            Inv = Matrix(inv(self.entries))
627        Inv.set("inverse", self)
628        # Set known properties
629        if "determinant" in self._computed:
630            Inv.set("determinant", 1.0 / self.determinant())
631        if "eigenvalues" in self._computed:
632            Inv.set("eigenvalues", np.flip(1.0 / self.eigenvalues()))
633        # Set the lin solver to the found inverse
634        self.set("lin_solver", lambda v: solve.invert(self, v))
635        return Inv
636
637    @_computable_property("lin_solver")
638    def lin_solver(self):
639        """Return a method which solves $Ax=b$ for $x$.
640        See `solve.automatic`."""
641        return solve.automatic(self, None, get_method=True)
642
643    @_computable_property("qr")
644    def qr(self):
645        """Compute the QR factorization of the matrix. See `factor.qr`."""
646        return factor.qr(self)
647
648    @_computable_property("svd")
649    def svd(self):
650        """Compute the SVD factorization of the matrix. See `factor.svd`."""
651        return factor.svd(self)
652
653    @_computable_property("transpose")
654    def transpose(self):
655        """Compute the transpose of the matrix."""
656        if self.is_sparse():
657            T = Matrix(self.entries.transpose())
658        else:
659            T = Matrix(self.entries.T)
660        # Copy properties that are identical for the transpose
661        T._computed = {prop: self._computed[prop] for prop in [
662            "complex", "cond_2", "determinant", "diagonal"
663            "eigenvalues", "eigenvectors", "invertible", "rank",
664            "sigmas", "symmetric"
665        ] if prop in self._computed}
666        T.set("transpose", self)
667        if "triangular_L" in self._computed:
668            T.set("triangular_U", self.is_triangular_L())
669        if "triangular_U" in self._computed:
670            T.set("triangular_L", self.is_triangular_U())
671        if "svd" in self._computed:
672            U, S, V = self.svd()
673            T.set("svd", (V.transpose(), S, U.transpose()))
674        return T
675
676    """Operation implementations."""
677
678    def __mul__(self, M):
679        """Multiplication by a constant."""
680        if type(M) is int or type(M) is float:
681            return Matrix(M * self.entries)
682        raise TypeError("Use of * to multiply by a {}.".format(type(M)))
683
684    def __rmul__(self, M):
685        return self.__mul__(M)
686
687    def __neg__(self):
688        return self.__mul__(-1)
689
690    def __matmul__(self, M):
691        """Matrix multiplication."""
692        if not isinstance(M, Matrix):
693            # Try to promote to Matrix object
694            M = Matrix(M)
695        if self.nc == M.nr:
696            return Matrix(self.entries @ M.entries)
697        raise MatrixSizeException(self, M, "multiply")
698
699    def __add__(self, M):
700        """Addition"""
701        if self.nc != M.nc or self.nr != M.nr:
702            raise MatrixSizeException(self, M, "addition")
703        if not isinstance(M, Matrix):
704            # Try to promote to Matrix object
705            M = Matrix(M)
706        return Matrix(self.entries + M.entries)
707
708    def __radd__(self, M):
709        return self.__add__(M)
710
711    def __sub__(self, M):
712        """Subtraction"""
713        if self.nc != M.nc or self.nr != M.nr:
714            raise MatrixSizeException(self, M, "subtraction")
715        if not isinstance(M, Matrix):
716            # Try to promote to Matrix object
717            M = Matrix(M)
718        return Matrix(self.entries - M.entries)
719
720    def __rsub__(self, M):
721        return (-self).__add__(M)
722
723
724class DiagonalMatrix(Matrix):
725    """A DiagonalMatrix is a Matrix whose only nonzero entries
726    are on the primary diagonal."""
727
728    def __init__(self, diagonal):
729        """Initialize a DiagonalMatrix from an array of the diagonal."""
730        if isinstance(diagonal, Vector):
731            diagonal = diagonal.entries
732        # Use the SciPy sparse diagonal implementation.
733        super().__init__(dia_array((diagonal, 0),
734                                   shape=(len(diagonal), len(diagonal))))
735        self.set("diagonal", True)
736
737    def __getitem__(self, number):
738        """Return the 1-indexed entry along the diagonal."""
739        if number < 0:
740            return self.entries.diagonal()[number]
741        return self.entries.diagonal()[number-1]
742
743    @_computable_property("invertible")
744    def is_invertible(self, tolerance=1e-10):
745        """`True` if there are no zeros on the diagonal."""
746        if "inverse" in self._computed:
747            return True
748        return np.isclose(
749            self.entries.diagonal(), 0, atol=tolerance).sum() == 0
750
751    @_computable_property("determinant")
752    def determinant(self):
753        """Compute the determinant as the product of the diagonal."""
754        return np.prod(self.entries.diagonal())
755
756    @_computable_property("eigenvalues")
757    def eigenvalues(self, tolerance=1e-10):
758        """Compute the eigenvalues, equivalent to sorting the diagonal."""
759        return np.sort(self.entries.diagonal())
760
761    @_computable_property("inverse")
762    def inverse(self):
763        """Compute the inverse as the reciprocal of the diagonal."""
764        if self.is_invertible():
765            return DiagonalMatrix(1.0 / self.entries.diagonal())
766        raise MatrixTypeException("inverse", "invertible")
767
768    @_computable_property("trace")
769    def trace(self):
770        """Compute the trace of the matrix."""
771        return np.sum(self.diagonal().entries)
772
773
774class Identity(DiagonalMatrix):
775    """The Identity matrix is a DiagonalMatrix of ones."""
776
777    def __init__(self, n):
778        """Initialize the Identity of size n."""
779        Matrix.__init__(self, identity(n))
780
781
782class Vector(Matrix):
783    """A Vector is a Matrix of size 1 in one dimension."""
784
785    def __init__(self, arr, row=False):
786        """Initialize a Vector object from a Numpy array.
787        If the array has only one dimension, defaults to being a column vector.
788        A row vector can be forced by setting row=True."""
789        if isinstance(arr, Matrix):
790            arr = arr.entries
791        arr = np.array(arr)
792        if len(arr.shape) == 1:
793            # One-dimensional array.
794            if row:
795                super().__init__(arr.reshape(1, len(arr)))
796            else:
797                super().__init__(arr.reshape(len(arr), 1))
798        else:
799            # 2-dimensional array, initialize the Matrix.
800            super().__init__(arr)
801
802        if min(self.nr, self.nc) != 1:
803            raise ValueError(
804                "A Vector cannot have size {}x{}".format(self.nr, self.nc))
805
806    def __len__(self):
807        """Return the number of entries in the Vector."""
808        return max(self.nr, self.nc)
809
810    def __getitem__(self, number):
811        """Returns the 1-indexed entry in the Vector."""
812        if number < 0:
813            return self.entries[number]
814        return self.entries[number-1]
815
816    def __repr__(self):
817        return "{}x{} Vector ({} properties)".format(
818            self.nr, self.nc, len(self._computed))
819
820    @_computable_property("unit")
821    def unit(self):
822        """Compute the unit vector in the direction of this Vector.
823        Computes the 2-norm."""
824        return Vector(self.entries / self.norm_2())