Repository: sympy/sympy

Issue:
QRdecomposition fails for some matrices
```
Hello,

on the current git version of sympy I got a strange failure of the
Matrix.QRdecomposition method:

===========================================================================
In [1]: A = Matrix([[1,1],[1,1.01]])

In [2]: A.QRdecomposition()
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)

/Users/voss/sympy/<ipython console> in <module>()

/Users/voss/sympy/sympy/matrices/matrices.py in QRdecomposition(self)
    966             R[j,j] = tmp.norm()
    967             Q[:,j] = tmp / R[j,j]
--> 968             assert Q[:,j].norm() == 1
    969             for i in range(j):
    970                 R[i,j] = Q[:,i].dot(self[:,j])

AssertionError: 
===========================================================================

This seems not to be a problem with the conditioning of A, since it works
in the following case (which has a worse conditioning):

===========================================================================
In [3]: A = Matrix([[1,1],[1,1.001]])

In [4]: A.QRdecomposition()
Out[4]: 
⎛⎡  ⎽⎽⎽                    ⎤                              ⎞
⎜⎢╲╱ 2                     ⎥  ⎡  ⎽⎽⎽              ⎽⎽⎽    ⎤⎟
⎜⎢─────  -0.707106781186548⎥, ⎢╲╱ 2      1.0005⋅╲╱ 2     ⎥⎟
⎜⎢  2                      ⎥  ⎢                          ⎥⎟
⎜⎢                         ⎥  ⎣  0    0.00070710678118647⎦⎟
⎜⎢  ⎽⎽⎽                    ⎥                              ⎟
⎜⎢╲╱ 2                     ⎥                              ⎟
⎜⎢─────  0.707106781186548 ⎥                              ⎟
⎝⎣  2                      ⎦                              ⎠
===========================================================================

I hope this helps,
Jochen
```

Original issue for #4520: http://code.google.com/p/sympy/issues/detail?id=1421
Original author: https://code.google.com/u/jochen.voss@gmail.com/
Original owner: https://code.google.com/u/jochen.voss@gmail.com/


Discussion:

Patch:
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 8ae120bb9d..b821c42845 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -7,7 +7,7 @@
 from sympy.utilities import sift
 from sympy.utilities.misc import filldedent
 
-from sympy.matrices.expressions.matexpr import MatrixExpr, ZeroMatrix, Identity, MatrixElement
+from sympy.matrices.expressions.matexpr import MatrixExpr, ZeroMatrix, Identity
 from sympy.matrices.expressions.matmul import MatMul
 from sympy.matrices.expressions.matadd import MatAdd
 from sympy.matrices.expressions.matpow import MatPow
@@ -234,24 +234,16 @@ def transpose(self):
 
     def _entry(self, i, j, **kwargs):
         # Find row entry
-        orig_i, orig_j = i, j
         for row_block, numrows in enumerate(self.rowblocksizes):
-            cmp = i < numrows
-            if cmp == True:
+            if (i < numrows) != False:
                 break
-            elif cmp == False:
+            else:
                 i -= numrows
-            elif row_block < self.blockshape[0] - 1:
-                # Can't tell which block and it's not the last one, return unevaluated
-                return MatrixElement(self, orig_i, orig_j)
         for col_block, numcols in enumerate(self.colblocksizes):
-            cmp = j < numcols
-            if cmp == True:
+            if (j < numcols) != False:
                 break
-            elif cmp == False:
+            else:
                 j -= numcols
-            elif col_block < self.blockshape[1] - 1:
-                return MatrixElement(self, orig_i, orig_j)
         return self.blocks[row_block, col_block][i, j]
 
     @property
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index c4dcd4aedc..2a1789fc9f 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -1036,16 +1036,11 @@ class OneMatrix(MatrixExpr):
     """
     Matrix whose all entries are ones.
     """
-    def __new__(cls, m, n, evaluate=False):
+    def __new__(cls, m, n):
         m, n = _sympify(m), _sympify(n)
         cls._check_dim(m)
         cls._check_dim(n)
 
-        if evaluate:
-            condition = Eq(m, 1) & Eq(n, 1)
-            if condition == True:
-                return Identity(1)
-
         obj = super(OneMatrix, cls).__new__(cls, m, n)
         return obj
 
@@ -1057,12 +1052,6 @@ def as_explicit(self):
         from sympy import ImmutableDenseMatrix
         return ImmutableDenseMatrix.ones(*self.shape)
 
-    def doit(self, **hints):
-        args = self.args
-        if hints.get('deep', True):
-            args = [a.doit(**hints) for a in args]
-        return self.func(*args, evaluate=True)
-
     def _eval_transpose(self):
         return OneMatrix(self.cols, self.rows)
 
diff --git a/sympy/matrices/expressions/matmul.py b/sympy/matrices/expressions/matmul.py
index d84f0b8365..8c9e6016b2 100644
--- a/sympy/matrices/expressions/matmul.py
+++ b/sympy/matrices/expressions/matmul.py
@@ -9,7 +9,7 @@
 
 from .inverse import Inverse
 from .matexpr import \
-    MatrixExpr, ShapeError, Identity, ZeroMatrix, OneMatrix, GenericIdentity
+    MatrixExpr, ShapeError, Identity, ZeroMatrix, GenericIdentity
 from .matpow import MatPow
 from .transpose import transpose
 from .permutation import PermutationMatrix
@@ -363,28 +363,8 @@ def combine_permutations(mul):
 
     return MatMul(*result)
 
-def combine_one_matrices(mul):
-    """
-    Combine products of OneMatrix
-
-    e.g. OneMatrix(2, 3) * OneMatrix(3, 4) -> 3 * OneMatrix(2, 4)
-    """
-    factor, args = mul.as_coeff_matrices()
-    new_args = [args[0]]
-
-    for B in args[1:]:
-        A = new_args[-1]
-        if not isinstance(A, OneMatrix) or not isinstance(B, OneMatrix):
-            new_args.append(B)
-            continue
-        new_args.pop()
-        new_args.append(OneMatrix(A.shape[0], B.shape[1]))
-        factor *= A.shape[1]
-
-    return newmul(factor, *new_args)
-
 rules = (
-    any_zeros, remove_ids, combine_one_matrices, combine_powers, unpack, rm_id(lambda x: x == 1),
+    any_zeros, remove_ids, combine_powers, unpack, rm_id(lambda x: x == 1),
     merge_explicit, factor_in_front, flatten, combine_permutations)
 
 canonicalize = exhaust(typed({MatMul: do_one(*rules)}))
diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index 0a94d81db3..34a35faed2 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -729,12 +729,6 @@ def size(self):
             return S.Infinity
         return Integer(abs(dif//self.step))
 
-    @property
-    def is_finite_set(self):
-        if self.start.is_integer and self.stop.is_integer:
-            return True
-        return self.size.is_finite
-
     def __nonzero__(self):
         return self.start != self.stop
 


Test patch:
diff --git a/sympy/matrices/expressions/tests/test_blockmatrix.py b/sympy/matrices/expressions/tests/test_blockmatrix.py
index 9bc99b0ca0..da17b9dbb5 100644
--- a/sympy/matrices/expressions/tests/test_blockmatrix.py
+++ b/sympy/matrices/expressions/tests/test_blockmatrix.py
@@ -192,6 +192,7 @@ def test_BlockDiagMatrix():
 def test_blockcut():
     A = MatrixSymbol('A', n, m)
     B = blockcut(A, (n/2, n/2), (m/2, m/2))
+    assert A[i, j] == B[i, j]
     assert B == BlockMatrix([[A[:n/2, :m/2], A[:n/2, m/2:]],
                              [A[n/2:, :m/2], A[n/2:, m/2:]]])
 
diff --git a/sympy/matrices/expressions/tests/test_indexing.py b/sympy/matrices/expressions/tests/test_indexing.py
index dc3f8dbae8..737786be25 100644
--- a/sympy/matrices/expressions/tests/test_indexing.py
+++ b/sympy/matrices/expressions/tests/test_indexing.py
@@ -1,7 +1,7 @@
 from sympy import (symbols, MatrixSymbol, MatPow, BlockMatrix, KroneckerDelta,
         Identity, ZeroMatrix, ImmutableMatrix, eye, Sum, Dummy, trace,
         Symbol)
-from sympy.testing.pytest import raises, XFAIL
+from sympy.testing.pytest import raises
 from sympy.matrices.expressions.matexpr import MatrixElement, MatrixExpr
 
 k, l, m, n = symbols('k l m n', integer=True)
@@ -83,72 +83,6 @@ def test_block_index():
     assert BI.as_explicit().equals(eye(6))
 
 
-def test_block_index_symbolic():
-    # Note that these matrices may be zero-sized and indices may be negative, which causes
-    # all naive simplifications given in the comments to be invalid
-    A1 = MatrixSymbol('A1', n, k)
-    A2 = MatrixSymbol('A2', n, l)
-    A3 = MatrixSymbol('A3', m, k)
-    A4 = MatrixSymbol('A4', m, l)
-    A = BlockMatrix([[A1, A2], [A3, A4]])
-    assert A[0, 0] == MatrixElement(A, 0, 0)  # Cannot be A1[0, 0]
-    assert A[n - 1, k - 1] == A1[n - 1, k - 1]
-    assert A[n, k] == A4[0, 0]
-    assert A[n + m - 1, 0] == MatrixElement(A, n + m - 1, 0)  # Cannot be A3[m - 1, 0]
-    assert A[0, k + l - 1] == MatrixElement(A, 0, k + l - 1)  # Cannot be A2[0, l - 1]
-    assert A[n + m - 1, k + l - 1] == MatrixElement(A, n + m - 1, k + l - 1)  # Cannot be A4[m - 1, l - 1]
-    assert A[i, j] == MatrixElement(A, i, j)
-    assert A[n + i, k + j] == MatrixElement(A, n + i, k + j)  # Cannot be A4[i, j]
-    assert A[n - i - 1, k - j - 1] == MatrixElement(A, n - i - 1, k - j - 1)  # Cannot be A1[n - i - 1, k - j - 1]
-
-
-def test_block_index_symbolic_nonzero():
-    # All invalid simplifications from test_block_index_symbolic() that become valid if all
-    # matrices have nonzero size and all indices are nonnegative
-    k, l, m, n = symbols('k l m n', integer=True, positive=True)
-    i, j = symbols('i j', integer=True, nonnegative=True)
-    A1 = MatrixSymbol('A1', n, k)
-    A2 = MatrixSymbol('A2', n, l)
-    A3 = MatrixSymbol('A3', m, k)
-    A4 = MatrixSymbol('A4', m, l)
-    A = BlockMatrix([[A1, A2], [A3, A4]])
-    assert A[0, 0] == A1[0, 0]
-    assert A[n + m - 1, 0] == A3[m - 1, 0]
-    assert A[0, k + l - 1] == A2[0, l - 1]
-    assert A[n + m - 1, k + l - 1] == A4[m - 1, l - 1]
-    assert A[i, j] == MatrixElement(A, i, j)
-    assert A[n + i, k + j] == A4[i, j]
-    assert A[n - i - 1, k - j - 1] == A1[n - i - 1, k - j - 1]
-    assert A[2 * n, 2 * k] == A4[n, k]
-
-
-def test_block_index_large():
-    n, m, k = symbols('n m k', integer=True, positive=True)
-    i = symbols('i', integer=True, nonnegative=True)
-    A1 = MatrixSymbol('A1', n, n)
-    A2 = MatrixSymbol('A2', n, m)
-    A3 = MatrixSymbol('A3', n, k)
-    A4 = MatrixSymbol('A4', m, n)
-    A5 = MatrixSymbol('A5', m, m)
-    A6 = MatrixSymbol('A6', m, k)
-    A7 = MatrixSymbol('A7', k, n)
-    A8 = MatrixSymbol('A8', k, m)
-    A9 = MatrixSymbol('A9', k, k)
-    A = BlockMatrix([[A1, A2, A3], [A4, A5, A6], [A7, A8, A9]])
-    assert A[n + i, n + i] == MatrixElement(A, n + i, n + i)
-
-
-@XFAIL
-def test_block_index_symbolic_fail():
-    # To make this work, symbolic matrix dimensions would need to be somehow assumed nonnegative
-    # even if the symbols aren't specified as such.  Then 2 * n < n would correctly evaluate to
-    # False in BlockMatrix._entry()
-    A1 = MatrixSymbol('A1', n, 1)
-    A2 = MatrixSymbol('A2', m, 1)
-    A = BlockMatrix([[A1], [A2]])
-    assert A[2 * n, 0] == A2[n, 0]
-
-
 def test_slicing():
     A.as_explicit()[0, :]  # does not raise an error
 
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index 3fa18bae64..1d2ec09441 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -62,7 +62,6 @@ def test_zero_matrix_creation():
 def test_one_matrix_creation():
     assert OneMatrix(2, 2)
     assert OneMatrix(0, 0)
-    assert Eq(OneMatrix(1, 1), Identity(1))
     raises(ValueError, lambda: OneMatrix(-1, 2))
     raises(ValueError, lambda: OneMatrix(2.0, 2))
     raises(ValueError, lambda: OneMatrix(2j, 2))
@@ -178,12 +177,6 @@ def test_OneMatrix_doit():
     assert isinstance(Unn.doit().rows, Mul)
 
 
-def test_OneMatrix_mul():
-    assert OneMatrix(n, m) * OneMatrix(m, k) == OneMatrix(n, k) * m
-    assert w * OneMatrix(1, 1) == w
-    assert OneMatrix(1, 1) * w.T == w.T
-
-
 def test_Identity():
     A = MatrixSymbol('A', n, m)
     i, j = symbols('i j')
diff --git a/sympy/matrices/tests/test_decompositions.py b/sympy/matrices/tests/test_decompositions.py
index adacd31d2b..8310c69f68 100644
--- a/sympy/matrices/tests/test_decompositions.py
+++ b/sympy/matrices/tests/test_decompositions.py
@@ -3,6 +3,8 @@
 from sympy.matrices import Matrix, zeros, eye, SparseMatrix
 from sympy.abc import x, y, z
 from sympy.testing.pytest import raises
+from sympy.testing.matrices import allclose
+
 
 def test_LUdecomp():
     testmat = Matrix([[0, 2, 5, 3],
@@ -215,6 +217,21 @@ def test_QR_trivial():
     assert R.is_upper
     assert A == Q*R
 
+
+def test_QR_float():
+    A = Matrix([[1, 1], [1, 1.01]])
+    Q, R = A.QRdecomposition()
+    assert allclose(Q * R, A)
+    assert allclose(Q * Q.T, Matrix.eye(2))
+    assert allclose(Q.T * Q, Matrix.eye(2))
+
+    A = Matrix([[1, 1], [1, 1.001]])
+    Q, R = A.QRdecomposition()
+    assert allclose(Q * R, A)
+    assert allclose(Q * Q.T, Matrix.eye(2))
+    assert allclose(Q.T * Q, Matrix.eye(2))
+
+
 def test_LUdecomposition_Simple_iszerofunc():
     # Test if callable passed to matrices.LUdecomposition_Simple() as iszerofunc keyword argument is used inside
     # matrices.LUdecomposition_Simple()
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index f89c571742..4d45df5175 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -446,28 +446,6 @@ def test_range_interval_intersection():
     assert Range(0).intersect(Interval(0.2, 0.8)) is S.EmptySet
     assert Range(0).intersect(Interval(-oo, oo)) is S.EmptySet
 
-def test_range_is_finite_set():
-    assert Range(-100, 100).is_finite_set is True
-    assert Range(2, oo).is_finite_set is False
-    assert Range(-oo, 50).is_finite_set is False
-    assert Range(-oo, oo).is_finite_set is False
-    assert Range(oo, -oo).is_finite_set is True
-    assert Range(0, 0).is_finite_set is True
-    assert Range(oo, oo).is_finite_set is True
-    assert Range(-oo, -oo).is_finite_set is True
-    n = Symbol('n', integer=True)
-    m = Symbol('m', integer=True)
-    assert Range(n, n + 49).is_finite_set is True
-    assert Range(n, 0).is_finite_set is True
-    assert Range(-3, n + 7).is_finite_set is True
-    assert Range(n, m).is_finite_set is True
-    assert Range(n + m, m - n).is_finite_set is True
-    assert Range(n, n + m + n).is_finite_set is True
-    assert Range(n, oo).is_finite_set is False
-    assert Range(-oo, n).is_finite_set is False
-    # assert Range(n, -oo).is_finite_set is True
-    # assert Range(oo, n).is_finite_set is True
-    # Above tests fail due to a (potential) bug in sympy.sets.fancysets.Range.size (See issue #18999)
 
 def test_Integers_eval_imageset():
     ans = ImageSet(Lambda(x, 2*x + Rational(3, 7)), S.Integers)
diff --git a/sympy/testing/matrices.py b/sympy/testing/matrices.py
new file mode 100644
index 0000000000..236a384366
--- /dev/null
+++ b/sympy/testing/matrices.py
@@ -0,0 +1,8 @@
+def allclose(A, B, rtol=1e-05, atol=1e-08):
+    if len(A) != len(B):
+        return False
+
+    for x, y in zip(A, B):
+        if abs(x-y) > atol + rtol * max(abs(x), abs(y)):
+            return False
+    return True