Repository: sympy/sympy
Pull request: Various block matrix improvements
#### References to other Issues or PRs
Fixes #2460

#### Brief description of what is fixed or changed
Various improvements regarding BlockMatrix:
- Improve inversion of 2x2 block matrices
- Add assumption handlers for OneMatrix
- Raise more correct and more specific exceptions for matrices
- Allow creating BlockDiagMatrix from mutable matrices
- Improve evaluation of BlockMatrix trace/determinant
- Fix BlockDiagMatrix with non-square blocks

#### Other comments
Please let me know if smaller pull requests are better.

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Problem statement:
Regression: can't make expressions with mutable BlockDiagMatrix's
This used to work in 0.7.3 :

``` Python
>>> from sympy import Matrix, BlockDiagMatrix
>>> from sympy.abc import x, y, z
>>> bdm1 = BlockDiagMatrix(Matrix([x]), Matrix([y]))
>>> bdm2 = BlockDiagMatrix(Matrix([y]), Matrix([z]))
>>> bdm1 + bdm2
TypeError: unhashable type: 'MutableDenseMatrix'
```

Hints from issue discussion:
The solution clearly is that BlockDiagMatrix should convert MutableMatrix to ImmutableMatrix.

I bisected to b085bab7427ce6b5c2ad7657f216445659b5c4a1

Commit message: Improve inversion of 2x2 block matrices

It previously tried to apply the formula even when A wasn't square.  It now
also considers invertibility of matrices to choose which formula to apply.

Patch to review:
diff --git a/sympy/assumptions/handlers/matrices.py b/sympy/assumptions/handlers/matrices.py
index 98762511fab7..e1f60893e7a8 100644
--- a/sympy/assumptions/handlers/matrices.py
+++ b/sympy/assumptions/handlers/matrices.py
@@ -78,6 +78,8 @@ def MatrixSymbol(expr, assumptions):
     def ZeroMatrix(expr, assumptions):
         return ask(Q.square(expr), assumptions)
 
+    OneMatrix = ZeroMatrix
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.symmetric(expr.arg), assumptions)
@@ -138,6 +140,10 @@ def MatrixSymbol(expr, assumptions):
 
     ZeroMatrix = staticmethod(CommonHandler.AlwaysFalse)
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.invertible(expr.arg), assumptions)
@@ -286,6 +292,10 @@ def MatPow(expr, assumptions):
 
     ZeroMatrix = staticmethod(CommonHandler.AlwaysFalse)
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.fullrank(expr.arg), assumptions)
@@ -337,6 +347,10 @@ def MatrixSymbol(expr, assumptions):
 
     ZeroMatrix = staticmethod(CommonHandler.AlwaysFalse)
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.positive_definite(expr.arg), assumptions)
@@ -386,6 +400,10 @@ def MatrixSymbol(expr, assumptions):
 
     Identity, ZeroMatrix = [staticmethod(CommonHandler.AlwaysTrue)]*2
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.lower_triangular(expr.arg), assumptions)
@@ -439,6 +457,10 @@ def MatrixSymbol(expr, assumptions):
 
     Identity, ZeroMatrix = [staticmethod(CommonHandler.AlwaysTrue)]*2
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.upper_triangular(expr.arg), assumptions)
@@ -502,6 +524,10 @@ def MatrixSymbol(expr, assumptions):
     def ZeroMatrix(expr, assumptions):
         return True
 
+    @staticmethod
+    def OneMatrix(expr, assumptions):
+        return expr.shape[0] == 1 and expr.shape[1] == 1
+
     @staticmethod
     def Transpose(expr, assumptions):
         return ask(Q.diagonal(expr.arg), assumptions)
@@ -567,7 +593,7 @@ def MatPow(expr, assumptions):
 
     HadamardProduct, Determinant, Trace, Transpose = [MatAdd]*4
 
-    ZeroMatrix, Identity = [staticmethod(CommonHandler.AlwaysTrue)]*2
+    ZeroMatrix, OneMatrix, Identity = [staticmethod(CommonHandler.AlwaysTrue)]*3
 
     MatMul = staticmethod(partial(MatMul_elements, Q.integer_elements,
                                                    Q.integer))
diff --git a/sympy/assumptions/tests/test_matrices.py b/sympy/assumptions/tests/test_matrices.py
index d5176e29c67e..0888c1e9a42e 100644
--- a/sympy/assumptions/tests/test_matrices.py
+++ b/sympy/assumptions/tests/test_matrices.py
@@ -1,6 +1,6 @@
 from sympy import Q, ask, Symbol, DiagMatrix, DiagonalMatrix
 from sympy.matrices.expressions import (MatrixSymbol, Identity, ZeroMatrix,
-        Trace, MatrixSlice, Determinant)
+        OneMatrix, Trace, MatrixSlice, Determinant)
 from sympy.matrices.expressions.factorizations import LofLU
 from sympy.testing.pytest import XFAIL
 
@@ -29,6 +29,8 @@ def test_invertible():
     assert ask(Q.invertible(X.I)) is True
     assert ask(Q.invertible(Identity(3))) is True
     assert ask(Q.invertible(ZeroMatrix(3, 3))) is False
+    assert ask(Q.invertible(OneMatrix(1, 1))) is True
+    assert ask(Q.invertible(OneMatrix(3, 3))) is False
     assert ask(Q.invertible(X), Q.fullrank(X) & Q.square(X))
 
 def test_singular():
@@ -58,6 +60,9 @@ def test_symmetric():
     assert ask(Q.symmetric(V1.T*(V1 + V2))) is True
     assert ask(Q.symmetric(V1.T*(V1 + V2) + A1x1)) is True
     assert ask(Q.symmetric(MatrixSlice(Y, (0, 1), (1, 2)))) is True
+    assert ask(Q.symmetric(Identity(3))) is True
+    assert ask(Q.symmetric(ZeroMatrix(3, 3))) is True
+    assert ask(Q.symmetric(OneMatrix(3, 3))) is True
 
 def _test_orthogonal_unitary(predicate):
     assert ask(predicate(X), predicate(X))
@@ -89,6 +94,8 @@ def test_fullrank():
     assert ask(Q.fullrank(X*Z), Q.fullrank(X) & Q.fullrank(Z)) is True
     assert ask(Q.fullrank(Identity(3))) is True
     assert ask(Q.fullrank(ZeroMatrix(3, 3))) is False
+    assert ask(Q.fullrank(OneMatrix(1, 1))) is True
+    assert ask(Q.fullrank(OneMatrix(3, 3))) is False
     assert ask(Q.invertible(X), ~Q.fullrank(X)) == False
 
 
@@ -107,6 +114,8 @@ def test_positive_definite():
     assert not ask(Q.positive_definite(Y.T*X*Y), Q.positive_definite(X))
     assert ask(Q.positive_definite(Identity(3))) is True
     assert ask(Q.positive_definite(ZeroMatrix(3, 3))) is False
+    assert ask(Q.positive_definite(OneMatrix(1, 1))) is True
+    assert ask(Q.positive_definite(OneMatrix(3, 3))) is False
     assert ask(Q.positive_definite(X + Z), Q.positive_definite(X) &
             Q.positive_definite(Z)) is True
     assert not ask(Q.positive_definite(-X), Q.positive_definite(X))
@@ -119,6 +128,11 @@ def test_triangular():
             Q.lower_triangular(Z)) is True
     assert ask(Q.lower_triangular(Identity(3))) is True
     assert ask(Q.lower_triangular(ZeroMatrix(3, 3))) is True
+    assert ask(Q.upper_triangular(ZeroMatrix(3, 3))) is True
+    assert ask(Q.lower_triangular(OneMatrix(1, 1))) is True
+    assert ask(Q.upper_triangular(OneMatrix(1, 1))) is True
+    assert ask(Q.lower_triangular(OneMatrix(3, 3))) is False
+    assert ask(Q.upper_triangular(OneMatrix(3, 3))) is False
     assert ask(Q.triangular(X), Q.unit_triangular(X))
     assert ask(Q.upper_triangular(X**3), Q.upper_triangular(X))
     assert ask(Q.lower_triangular(X**3), Q.lower_triangular(X))
@@ -128,6 +142,8 @@ def test_diagonal():
     assert ask(Q.diagonal(X + Z.T + Identity(2)), Q.diagonal(X) &
                Q.diagonal(Z)) is True
     assert ask(Q.diagonal(ZeroMatrix(3, 3)))
+    assert ask(Q.diagonal(OneMatrix(1, 1))) is True
+    assert ask(Q.diagonal(OneMatrix(3, 3))) is False
     assert ask(Q.lower_triangular(X) & Q.upper_triangular(X), Q.diagonal(X))
     assert ask(Q.diagonal(X), Q.lower_triangular(X) & Q.upper_triangular(X))
     assert ask(Q.symmetric(X), Q.diagonal(X))
@@ -214,6 +230,7 @@ def test_matrix_element_sets():
     assert ask(Q.complex(X[1, 2]), Q.complex_elements(X))
     assert ask(Q.integer_elements(Identity(3)))
     assert ask(Q.integer_elements(ZeroMatrix(3, 3)))
+    assert ask(Q.integer_elements(OneMatrix(3, 3)))
     from sympy.matrices.expressions.fourier import DFT
     assert ask(Q.complex_elements(DFT(3)))
 
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 0f55936effdd..a8fda1eaef7a 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -1,7 +1,8 @@
 from __future__ import print_function, division
 
 from sympy import ask, Q
-from sympy.core import Basic, Add
+from sympy.core import Basic, Add, Mul, S
+from sympy.matrices.common import NonInvertibleMatrixError
 from sympy.strategies import typed, exhaust, condition, do_one, unpack
 from sympy.strategies.traverse import bottom_up
 from sympy.utilities import sift
@@ -12,11 +13,11 @@
 from sympy.matrices.expressions.matadd import MatAdd
 from sympy.matrices.expressions.matpow import MatPow
 from sympy.matrices.expressions.transpose import Transpose, transpose
-from sympy.matrices.expressions.trace import Trace
+from sympy.matrices.expressions.trace import trace
 from sympy.matrices.expressions.determinant import det, Determinant
 from sympy.matrices.expressions.slice import MatrixSlice
 from sympy.matrices.expressions.inverse import Inverse
-from sympy.matrices import Matrix, ShapeError
+from sympy.matrices import Matrix, ShapeError, MutableMatrix
 from sympy.functions.elementary.complexes import re, im
 
 class BlockMatrix(MatrixExpr):
@@ -186,12 +187,14 @@ def _eval_transpose(self):
 
     def _eval_trace(self):
         if self.rowblocksizes == self.colblocksizes:
-            return Add(*[Trace(self.blocks[i, i])
+            return Add(*[trace(self.blocks[i, i])
                         for i in range(self.blockshape[0])])
         raise NotImplementedError(
             "Can't perform trace of irregular blockshape")
 
     def _eval_determinant(self):
+        if self.blockshape == (1, 1):
+            return det(self.blocks[0, 0])
         if self.blockshape == (2, 2):
             [[A, B],
              [C, D]] = self.blocks.tolist()
@@ -306,6 +309,7 @@ class BlockDiagMatrix(BlockMatrix):
     sympy.matrices.dense.diag
     """
     def __new__(cls, *mats):
+        mats = [m.as_immutable() if isinstance(m, MutableMatrix) else m for m in mats]
         return Basic.__new__(BlockDiagMatrix, *mats)
 
     @property
@@ -339,8 +343,22 @@ def rowblocksizes(self):
     def colblocksizes(self):
         return [block.cols for block in self.args]
 
+    def _all_square_blocks(self):
+        """Returns true if all blocks are square"""
+        return all(mat.is_square for mat in self.args)
+
+    def _eval_determinant(self):
+        if self._all_square_blocks():
+            return Mul(*[det(mat) for mat in self.args])
+        # At least one block is non-square.  Since the entire matrix must be square we know there must
+        # be at least two blocks in this matrix, in which case the entire matrix is necessarily rank-deficient
+        return S.Zero
+
     def _eval_inverse(self, expand='ignored'):
-        return BlockDiagMatrix(*[mat.inverse() for mat in self.args])
+        if self._all_square_blocks():
+            return BlockDiagMatrix(*[mat.inverse() for mat in self.args])
+        # See comment in _eval_determinant()
+        raise NonInvertibleMatrixError('Matrix det == 0; not invertible.')
 
     def _eval_transpose(self):
         return BlockDiagMatrix(*[mat.transpose() for mat in self.args])
@@ -526,7 +544,7 @@ def bc_transpose(expr):
 
 def bc_inverse(expr):
     if isinstance(expr.arg, BlockDiagMatrix):
-        return expr._eval_inverse()
+        return expr.inverse()
 
     expr2 = blockinverse_1x1(expr)
     if expr != expr2:
@@ -540,15 +558,36 @@ def blockinverse_1x1(expr):
     return expr
 
 def blockinverse_2x2(expr):
-    if isinstance(expr.arg, BlockMatrix) and expr.arg.blockshape == (2, 2):
+    if isinstance(expr.arg, BlockMatrix) and expr.arg.blockshape == (2, 2) and expr.arg.blocks[0, 0].is_square:
         # Cite: The Matrix Cookbook Section 9.1.3
         [[A, B],
          [C, D]] = expr.arg.blocks.tolist()
 
-        return BlockMatrix([[ (A - B*D.I*C).I,  (-A).I*B*(D - C*A.I*B).I],
-                            [-(D - C*A.I*B).I*C*A.I,     (D - C*A.I*B).I]])
-    else:
-        return expr
+        # Use one or the other formula, depending on whether A or D is known to be invertible or at least not known
+        # to not be invertible.  Note that invertAbility of the other expressions M is not checked.
+        A_invertible = ask(Q.invertible(A))
+        D_invertible = ask(Q.invertible(D))
+        if A_invertible == True:
+            invert_A = True
+        elif D_invertible == True:
+            invert_A = False
+        elif A_invertible != False:
+            invert_A = True
+        elif D_invertible != False:
+            invert_A = False
+        else:
+            invert_A = True
+
+        if invert_A:
+            AI = A.I
+            MI = (D - C * AI * B).I
+            return BlockMatrix([[AI + AI * B * MI * C * AI, -AI * B * MI], [-MI * C * AI, MI]])
+        else:
+            DI = D.I
+            MI = (A - B * DI * C).I
+            return BlockMatrix([[MI, -MI * B * DI], [-DI * C * MI, DI + DI * C * MI * B * DI]])
+
+    return expr
 
 def deblock(B):
     """ Flatten a BlockMatrix of BlockMatrices """
diff --git a/sympy/matrices/expressions/determinant.py b/sympy/matrices/expressions/determinant.py
index f0b568510030..9df21dbcf86b 100644
--- a/sympy/matrices/expressions/determinant.py
+++ b/sympy/matrices/expressions/determinant.py
@@ -1,7 +1,7 @@
 from __future__ import print_function, division
 
 from sympy import Basic, Expr, S, sympify
-from .matexpr import ShapeError
+from sympy.matrices.common import NonSquareMatrixError
 
 
 class Determinant(Expr):
@@ -26,7 +26,7 @@ def __new__(cls, mat):
             raise TypeError("Input to Determinant, %s, not a matrix" % str(mat))
 
         if not mat.is_square:
-            raise ShapeError("Det of a non-square matrix")
+            raise NonSquareMatrixError("Det of a non-square matrix")
 
         return Basic.__new__(cls, mat)
 
diff --git a/sympy/matrices/expressions/inverse.py b/sympy/matrices/expressions/inverse.py
index 7bfe41299220..97d7c9f7a4a6 100644
--- a/sympy/matrices/expressions/inverse.py
+++ b/sympy/matrices/expressions/inverse.py
@@ -3,7 +3,7 @@
 from sympy.core.sympify import _sympify
 from sympy.core import S, Basic
 
-from sympy.matrices.expressions.matexpr import ShapeError
+from sympy.matrices.common import NonSquareMatrixError
 from sympy.matrices.expressions.matpow import MatPow
 
 
@@ -41,7 +41,7 @@ def __new__(cls, mat, exp=S.NegativeOne):
         if not mat.is_Matrix:
             raise TypeError("mat should be a matrix")
         if not mat.is_square:
-            raise ShapeError("Inverse of non-square matrix %s" % mat)
+            raise NonSquareMatrixError("Inverse of non-square matrix %s" % mat)
         return Basic.__new__(cls, mat, exp)
 
     @property
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index c4dcd4aedcf8..155f72d8b225 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -12,7 +12,7 @@
 from sympy.core.sympify import SympifyError, _sympify
 from sympy.functions import conjugate, adjoint
 from sympy.functions.special.tensor_functions import KroneckerDelta
-from sympy.matrices import ShapeError
+from sympy.matrices.common import ShapeError, NonSquareMatrixError, NonInvertibleMatrixError
 from sympy.simplify import simplify
 from sympy.utilities.misc import filldedent
 
@@ -129,7 +129,7 @@ def __rmatmul__(self, other):
     @call_highest_priority('__rpow__')
     def __pow__(self, other):
         if not self.is_square:
-            raise ShapeError("Power of non-square matrix %s" % self)
+            raise NonSquareMatrixError("Power of non-square matrix %s" % self)
         elif self.is_Identity:
             return self
         elif other == S.Zero:
@@ -265,6 +265,8 @@ def T(self):
         return self.transpose()
 
     def inverse(self):
+        if not self.is_square:
+            raise NonSquareMatrixError('Inverse of non-square matrix')
         return self._eval_inverse()
 
     def inv(self):
@@ -969,11 +971,11 @@ def shape(self):
     @call_highest_priority('__rpow__')
     def __pow__(self, other):
         if other != 1 and not self.is_square:
-            raise ShapeError("Power of non-square matrix %s" % self)
+            raise NonSquareMatrixError("Power of non-square matrix %s" % self)
         if other == 0:
             return Identity(self.rows)
         if other < 1:
-            raise ValueError("Matrix det == 0; not invertible.")
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible")
         return self
 
     def _eval_transpose(self):
@@ -985,6 +987,9 @@ def _eval_trace(self):
     def _eval_determinant(self):
         return S.Zero
 
+    def _eval_inverse(self):
+        raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+
     def conjugate(self):
         return self
 
@@ -1069,8 +1074,13 @@ def _eval_transpose(self):
     def _eval_trace(self):
         return S.One*self.rows
 
+    def _is_1x1(self):
+        """Returns true if the matrix is known to be 1x1"""
+        shape = self.shape
+        return Eq(shape[0], 1) & Eq(shape[1], 1)
+
     def _eval_determinant(self):
-        condition = Eq(self.shape[0], 1) & Eq(self.shape[1], 1)
+        condition = self._is_1x1()
         if condition == True:
             return S.One
         elif condition == False:
@@ -1079,6 +1089,15 @@ def _eval_determinant(self):
             from sympy import Determinant
             return Determinant(self)
 
+    def _eval_inverse(self):
+        condition = self._is_1x1()
+        if condition == True:
+            return Identity(1)
+        elif condition == False:
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+        else:
+            return Inverse(self)
+
     def conjugate(self):
         return self
 
diff --git a/sympy/matrices/expressions/matpow.py b/sympy/matrices/expressions/matpow.py
index ca372a38ba78..d64746d0fd71 100644
--- a/sympy/matrices/expressions/matpow.py
+++ b/sympy/matrices/expressions/matpow.py
@@ -1,9 +1,11 @@
 from __future__ import print_function, division
 
-from .matexpr import MatrixExpr, ShapeError, Identity, ZeroMatrix
+from sympy.matrices.common import NonSquareMatrixError
+from .matexpr import MatrixExpr, Identity, ZeroMatrix
 from sympy.core import S
 from sympy.core.sympify import _sympify
 from sympy.matrices import MatrixBase
+from sympy.matrices.common import NonInvertibleMatrixError
 
 from .permutation import PermutationMatrix
 
@@ -40,7 +42,7 @@ def _entry(self, i, j, **kwargs):
         if isinstance(A, MatPow):
             # We still have a MatPow, make an explicit MatMul out of it.
             if not A.base.is_square:
-                raise ShapeError("Power of non-square matrix %s" % A.base)
+                raise NonSquareMatrixError("Power of non-square matrix %s" % A.base)
             elif A.exp.is_Integer and A.exp.is_positive:
                 A = MatMul(*[A.base for k in range(A.exp)])
             #elif A.exp.is_Integer and self.exp.is_negative:
@@ -74,7 +76,7 @@ def doit(self, **kwargs):
                 return base.func(Identity(base.shape[0]))
             return Identity(base.shape[0])
         elif isinstance(base, ZeroMatrix) and exp.is_negative:
-            raise ValueError("Matrix determinant is 0, not invertible.")
+            raise NonInvertibleMatrixError("Matrix determinant is 0, not invertible")
         elif isinstance(base, (Identity, ZeroMatrix)):
             return base
         elif isinstance(base, PermutationMatrix):
diff --git a/sympy/matrices/expressions/tests/test_blockmatrix.py b/sympy/matrices/expressions/tests/test_blockmatrix.py
index 07d8d3ce0e5b..d92fbe5a36a6 100644
--- a/sympy/matrices/expressions/tests/test_blockmatrix.py
+++ b/sympy/matrices/expressions/tests/test_blockmatrix.py
@@ -1,10 +1,12 @@
+from sympy import Trace
 from sympy.testing.pytest import raises
 from sympy.matrices.expressions.blockmatrix import (
     block_collapse, bc_matmul, bc_block_plus_ident, BlockDiagMatrix,
     BlockMatrix, bc_dist, bc_matadd, bc_transpose, bc_inverse,
     blockcut, reblock_2x2, deblock)
 from sympy.matrices.expressions import (MatrixSymbol, Identity,
-        Inverse, trace, Transpose, det, ZeroMatrix)
+        Inverse, trace, Transpose, det, ZeroMatrix, OneMatrix)
+from sympy.matrices.common import NonInvertibleMatrixError
 from sympy.matrices import (
     Matrix, ImmutableMatrix, ImmutableSparseMatrix)
 from sympy.core import Tuple, symbols, Expr
@@ -119,6 +121,7 @@ def test_BlockMatrix_trace():
     A, B, C, D = [MatrixSymbol(s, 3, 3) for s in 'ABCD']
     X = BlockMatrix([[A, B], [C, D]])
     assert trace(X) == trace(A) + trace(D)
+    assert trace(BlockMatrix([ZeroMatrix(n, n)])) == 0
 
 def test_BlockMatrix_Determinant():
     A, B, C, D = [MatrixSymbol(s, 3, 3) for s in 'ABCD']
@@ -128,6 +131,8 @@ def test_BlockMatrix_Determinant():
         assert det(X) == det(A) * det(D - C*A.I*B)
 
     assert isinstance(det(X), Expr)
+    assert det(BlockMatrix([A])) == det(A)
+    assert det(BlockMatrix([ZeroMatrix(n, n)])) == 0
 
 def test_squareBlockMatrix():
     A = MatrixSymbol('A', n, n)
@@ -147,9 +152,6 @@ def test_squareBlockMatrix():
     assert (X * MatrixSymbol('Q', n + m, n + m)).is_MatMul
 
     assert block_collapse(Y.I) == A.I
-    assert block_collapse(X.inverse()) == BlockMatrix([
-        [(-B*D.I*C + A).I, -A.I*B*(D + -C*A.I*B).I],
-        [-(D - C*A.I*B).I*C*A.I, (D - C*A.I*B).I]])
 
     assert isinstance(X.inverse(), Inverse)
 
@@ -158,6 +160,34 @@ def test_squareBlockMatrix():
     Z = BlockMatrix([[Identity(n), B], [C, D]])
     assert not Z.is_Identity
 
+def test_BlockMatrix_inverse():
+    A = MatrixSymbol('A', n, m)
+    B = MatrixSymbol('B', n, n)
+    C = MatrixSymbol('C', m, m)
+    D = MatrixSymbol('D', m, n)
+    X = BlockMatrix([[A, B], [C, D]])
+    assert X.is_square
+    assert isinstance(block_collapse(X.inverse()), Inverse)  # Can't inverse when A, D aren't square
+
+    # test code path for non-invertible D matrix
+    A = MatrixSymbol('A', n, n)
+    B = MatrixSymbol('B', n, m)
+    C = MatrixSymbol('C', m, n)
+    D = OneMatrix(m, m)
+    X = BlockMatrix([[A, B], [C, D]])
+    assert block_collapse(X.inverse()) == BlockMatrix([
+        [A.I + A.I * B * (D - C * A.I * B).I * C * A.I, -A.I * B * (D - C * A.I * B).I],
+        [-(D - C * A.I * B).I * C * A.I, (D - C * A.I * B).I],
+    ])
+
+    # test code path for non-invertible A matrix
+    A = OneMatrix(n, n)
+    D = MatrixSymbol('D', m, m)
+    X = BlockMatrix([[A, B], [C, D]])
+    assert block_collapse(X.inverse()) == BlockMatrix([
+        [(A - B * D.I * C).I, -(A - B * D.I * C).I * B * D.I],
+        [-D.I * C * (A - B * D.I * C).I, D.I + D.I * C * (A - B * D.I * C).I * B * D.I],
+    ])
 
 def test_BlockDiagMatrix():
     A = MatrixSymbol('A', n, n)
@@ -191,6 +221,58 @@ def test_BlockDiagMatrix():
     assert (X._blockmul(M)).is_MatMul
     assert (X._blockadd(M)).is_MatAdd
 
+def test_BlockDiagMatrix_nonsquare():
+    A = MatrixSymbol('A', n, m)
+    B = MatrixSymbol('B', k, l)
+    X = BlockDiagMatrix(A, B)
+    assert X.shape == (n + k, m + l)
+    assert X.shape == (n + k, m + l)
+    assert X.rowblocksizes == [n, k]
+    assert X.colblocksizes == [m, l]
+    C = MatrixSymbol('C', n, m)
+    D = MatrixSymbol('D', k, l)
+    Y = BlockDiagMatrix(C, D)
+    assert block_collapse(X + Y) == BlockDiagMatrix(A + C, B + D)
+    assert block_collapse(X * Y.T) == BlockDiagMatrix(A * C.T, B * D.T)
+    raises(NonInvertibleMatrixError, lambda: BlockDiagMatrix(A, C.T).inverse())
+
+def test_BlockDiagMatrix_determinant():
+    A = MatrixSymbol('A', n, n)
+    B = MatrixSymbol('B', m, m)
+    assert det(BlockDiagMatrix()) == 1
+    assert det(BlockDiagMatrix(A)) == det(A)
+    assert det(BlockDiagMatrix(A, B)) == det(A) * det(B)
+
+    # non-square blocks
+    C = MatrixSymbol('C', m, n)
+    D = MatrixSymbol('D', n, m)
+    assert det(BlockDiagMatrix(C, D)) == 0
+
+def test_BlockDiagMatrix_trace():
+    assert trace(BlockDiagMatrix()) == 0
+    assert trace(BlockDiagMatrix(ZeroMatrix(n, n))) == 0
+    A = MatrixSymbol('A', n, n)
+    assert trace(BlockDiagMatrix(A)) == trace(A)
+    B = MatrixSymbol('B', m, m)
+    assert trace(BlockDiagMatrix(A, B)) == trace(A) + trace(B)
+
+    # non-square blocks
+    C = MatrixSymbol('C', m, n)
+    D = MatrixSymbol('D', n, m)
+    assert isinstance(trace(BlockDiagMatrix(C, D)), Trace)
+
+def test_BlockDiagMatrix_transpose():
+    A = MatrixSymbol('A', n, m)
+    B = MatrixSymbol('B', k, l)
+    assert transpose(BlockDiagMatrix()) == BlockDiagMatrix()
+    assert transpose(BlockDiagMatrix(A)) == BlockDiagMatrix(A.T)
+    assert transpose(BlockDiagMatrix(A, B)) == BlockDiagMatrix(A.T, B.T)
+
+def test_issue_2460():
+    bdm1 = BlockDiagMatrix(Matrix([i]), Matrix([j]))
+    bdm2 = BlockDiagMatrix(Matrix([k]), Matrix([l]))
+    assert block_collapse(bdm1 + bdm2) == BlockDiagMatrix(Matrix([i + k]), Matrix([j + l]))
+
 def test_blockcut():
     A = MatrixSymbol('A', n, m)
     B = blockcut(A, (n/2, n/2), (m/2, m/2))
diff --git a/sympy/matrices/expressions/tests/test_inverse.py b/sympy/matrices/expressions/tests/test_inverse.py
index 2fad7f110614..4fe753870f7f 100644
--- a/sympy/matrices/expressions/tests/test_inverse.py
+++ b/sympy/matrices/expressions/tests/test_inverse.py
@@ -1,6 +1,7 @@
 from sympy.core import symbols, S
-from sympy.matrices.expressions import MatrixSymbol, Inverse, MatPow
-from sympy.matrices import eye, Identity, ShapeError
+from sympy.matrices.expressions import MatrixSymbol, Inverse, MatPow, ZeroMatrix, OneMatrix
+from sympy.matrices.common import NonSquareMatrixError, NonInvertibleMatrixError
+from sympy.matrices import eye, Identity
 from sympy.testing.pytest import raises
 from sympy import refine, Q
 
@@ -13,9 +14,6 @@
 
 
 def test_inverse():
-    raises(ShapeError, lambda: Inverse(A))
-    raises(ShapeError, lambda: Inverse(A*B))
-
     assert Inverse(C).args == (C, S.NegativeOne)
     assert Inverse(C).shape == (n, n)
     assert Inverse(A*E).shape == (n, n)
@@ -41,6 +39,16 @@ def test_inverse():
     assert Inverse(eye(3)).doit() == eye(3)
     assert Inverse(eye(3)).doit(deep=False) == eye(3)
 
+    assert OneMatrix(1, 1).I == Identity(1)
+    assert isinstance(OneMatrix(n, n).I, Inverse)
+
+def test_inverse_non_invertible():
+    raises(NonSquareMatrixError, lambda: Inverse(A))
+    raises(NonSquareMatrixError, lambda: Inverse(A*B))
+    raises(NonSquareMatrixError, lambda: ZeroMatrix(n, m).I)
+    raises(NonInvertibleMatrixError, lambda: ZeroMatrix(n, n).I)
+    raises(NonSquareMatrixError, lambda: OneMatrix(n, m).I)
+    raises(NonInvertibleMatrixError, lambda: OneMatrix(2, 2).I)
 
 def test_refine():
     assert refine(C.I, Q.orthogonal(C)) == C.T
diff --git a/sympy/matrices/expressions/trace.py b/sympy/matrices/expressions/trace.py
index 3beccefd35a7..5fa6985badf7 100644
--- a/sympy/matrices/expressions/trace.py
+++ b/sympy/matrices/expressions/trace.py
@@ -2,7 +2,7 @@
 
 from sympy import Basic, Expr, sympify, S
 from sympy.matrices.matrices import MatrixBase
-from .matexpr import ShapeError
+from sympy.matrices.common import NonSquareMatrixError
 
 
 class Trace(Expr):
@@ -28,7 +28,7 @@ def __new__(cls, mat):
             raise TypeError("input to Trace, %s, is not a matrix" % str(mat))
 
         if not mat.is_square:
-            raise ShapeError("Trace of a non-square matrix")
+            raise NonSquareMatrixError("Trace of a non-square matrix")
 
         return Basic.__new__(cls, mat)
 
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 69be53842860..cfc55bf9fde1 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -2141,7 +2141,7 @@ def vech(self, diagonal=True, check_symmetry=True):
 
         c = self.cols
         if c != self.rows:
-            raise ShapeError("Matrix must be square")
+            raise NonSquareMatrixError("Matrix must be square")
         if check_symmetry:
             self.simplify()
             if self != self.transpose():
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index ecc408c22fd9..00ab5320e9ed 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -41,6 +41,7 @@
     separatevars)
 from sympy.simplify.sqrtdenest import sqrt_depth
 from sympy.simplify.fu import TR1
+from sympy.matrices.common import NonInvertibleMatrixError
 from sympy.matrices import Matrix, zeros
 from sympy.polys import roots, cancel, factor, Poly, degree
 from sympy.polys.polyerrors import GeneratorsNeeded, PolynomialError
@@ -2645,7 +2646,7 @@ def inv_quick(M):
     n = M.rows
     d = det(M)
     if d == S.Zero:
-        raise ValueError("Matrix det == 0; not invertible.")
+        raise NonInvertibleMatrixError("Matrix det == 0; not invertible")
     ret = zeros(n)
     s1 = -1
     for i in range(n):
diff --git a/sympy/stats/stochastic_process_types.py b/sympy/stats/stochastic_process_types.py
index 5d6b5e70a335..e29fc349be2e 100644
--- a/sympy/stats/stochastic_process_types.py
+++ b/sympy/stats/stochastic_process_types.py
@@ -5,7 +5,7 @@
                    Lambda, Mul, Dummy, IndexedBase, Add,
                    linsolve, eye, Or, Not, Intersection,
                    Union, Expr, Function, exp, cacheit,
-                   Ge, Piecewise, Symbol)
+                   Ge, Piecewise, Symbol, NonSquareMatrixError)
 from sympy.core.relational import Relational
 from sympy.logic.boolalg import Boolean
 from sympy.stats.joint_rv import JointDistributionHandmade, JointDistribution
@@ -80,7 +80,7 @@ def _matrix_checks(matrix):
         raise TypeError("Transition probabilities either should "
                             "be a Matrix or a MatrixSymbol.")
     if matrix.shape[0] != matrix.shape[1]:
-        raise ValueError("%s is not a square matrix"%(matrix))
+        raise NonSquareMatrixError("%s is not a square matrix"%(matrix))
     if isinstance(matrix, Matrix):
         matrix = ImmutableMatrix(matrix.tolist())
     return matrix


Code review comments:
diff_hunk: @@ -306,6 +309,7 @@ class BlockDiagMatrix(BlockMatrix):
     sympy.matrices.dense.diag
     """
     def __new__(cls, *mats):
+        mats = [m.as_immutable() if isinstance(m, MutableMatrix) else m for m in mats]
path: sympy/matrices/expressions/blockmatrix.py
text: @user1:
It's better to use _sympify here.