Repository: sympy/sympy

Issue:
gram schmidt bug `def _orthogonalize(cls, *vecs, **kwargs):` at line 107 + fix
I tried to orthogonalise 2 complex vectors but it gave me the wrong answer but I found out where the mistake was:
it used to say

- `def project(a, b):
        return b * (a.dot(b, hermitian=True) / b.dot(b, hermitian=True))`
But it should be
- `def project(a, b):
        return b * (b.dot(a, hermitian=True) / b.dot(b, hermitian=True))`

(In the file subspaces.py in the function: `def _orthogonalize(cls, *vecs, **kwargs):` at line 107.)

I found the mistake when orthogonalising 
`v1 = Matrix([1,I,0])
v2 = Matrix([I,0,-I])`
with the code in the first bullet point you get this:

`v1 = Matrix([1,I,0])
v2 = Matrix([I,0,-I])
W = [v1,v2]
GS = GramSchmidt(W)`

with GS= [Matrix([
[1],
[I],
[0]]), Matrix([
[3*I/2],
[ -1/2],
[   -I]])]
but these two vectors multiplied with the standard complex inner product you get 2i which is unequal to zero.

And with code in the second bullet point you do get zero hence they are orthogonal. 
In the slides of my university textbook the formulas correspond to code in the second bullet point:
![img](https://user-images.githubusercontent.com/40834779/97086834-16156000-1626-11eb-85df-0de8297e0f2e.png)
The formula I denoted in red being the formula I'm refering to. I'm a little unfamiliar with github, but I do want to contribute by making this software better since I really like sympy.

I hope what I wrote down is sort of readable and convincing enough for someone to fix this problem.


Discussion:

Patch:
diff --git a/sympy/matrices/decompositions.py b/sympy/matrices/decompositions.py
index 0d59313c3b..5330dab7c3 100644
--- a/sympy/matrices/decompositions.py
+++ b/sympy/matrices/decompositions.py
@@ -1123,6 +1123,45 @@ def _LUdecompositionFF(M):
     return P, L, DD, U
 
 
+def _QRdecomposition_optional(M, normalize=True):
+    def dot(u, v):
+        return u.dot(v, hermitian=True)
+
+    dps = _get_intermediate_simp(expand_mul, expand_mul)
+
+    A = M.as_mutable()
+    ranked = list()
+
+    Q = A
+    R = A.zeros(A.cols)
+
+    for j in range(A.cols):
+        for i in range(j):
+            if Q[:, i].is_zero_matrix:
+                continue
+
+            R[i, j] = dot(Q[:, i], Q[:, j]) / dot(Q[:, i], Q[:, i])
+            R[i, j] = dps(R[i, j])
+            Q[:, j] -= Q[:, i] * R[i, j]
+
+        Q[:, j] = dps(Q[:, j])
+        if Q[:, j].is_zero_matrix is False:
+            ranked.append(j)
+            R[j, j] = M.one
+
+    Q = Q.extract(range(Q.rows), ranked)
+    R = R.extract(ranked, range(R.cols))
+
+    if normalize:
+        # Normalization
+        for i in range(Q.cols):
+            norm = Q[:, i].norm()
+            Q[:, i] /= norm
+            R[i, :] *= norm
+
+    return M.__class__(Q), M.__class__(R)
+
+
 def _QRdecomposition(M):
     r"""Returns a QR decomposition.
 
@@ -1303,41 +1342,4 @@ def _QRdecomposition(M):
     LUdecomposition
     QRsolve
     """
-
-    dps    = _get_intermediate_simp(expand_mul, expand_mul)
-    mat    = M.as_mutable()
-    n      = mat.rows
-    m      = mat.cols
-    ranked = list()
-
-    # Pad with additional rows to make wide matrices square
-    # nOrig keeps track of original size so zeros can be trimmed from Q
-    if n < m:
-        nOrig = n
-        n     = m
-        mat   = mat.col_join(mat.zeros(n - nOrig, m))
-    else:
-        nOrig = n
-
-    Q, R = mat.zeros(n, m), mat.zeros(m)
-
-    for j in range(m):  # for each column vector
-        tmp = mat[:, j]  # take original v
-
-        for i in range(j):
-            # subtract the project of mat on new vector
-            R[i, j]  = dps(Q[:, i].dot(mat[:, j], hermitian=True))
-            tmp     -= Q[:, i] * R[i, j]
-
-        tmp = dps(tmp)
-
-        # normalize it
-        R[j, j] = tmp.norm()
-
-        if not R[j, j].is_zero:
-            ranked.append(j)
-            Q[:, j] = tmp / R[j, j]
-
-    Q = Q.extract(range(nOrig), ranked)
-    R = R.extract(ranked, range(R.cols))
-    return M.__class__(Q), M.__class__(R)
+    return _QRdecomposition_optional(M, normalize=True)
diff --git a/sympy/matrices/subspaces.py b/sympy/matrices/subspaces.py
index b5d039940f..846c510bc0 100644
--- a/sympy/matrices/subspaces.py
+++ b/sympy/matrices/subspaces.py
@@ -1,5 +1,3 @@
-from sympy.core.compatibility import reduce
-
 from .utilities import _iszero
 
 
@@ -152,39 +150,25 @@ def _orthogonalize(cls, *vecs, normalize=False, rankcheck=False):
 
     .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
     """
+    from .decompositions import _QRdecomposition_optional
 
-    def project(a, b):
-        return b * (a.dot(b, hermitian=True) / b.dot(b, hermitian=True))
-
-    def perp_to_subspace(vec, basis):
-        """projects vec onto the subspace given
-        by the orthogonal basis ``basis``"""
+    if not vecs:
+        return []
 
-        components = [project(vec, b) for b in basis]
+    all_row_vecs = (vecs[0].rows == 1)
 
-        if len(basis) == 0:
-            return vec
+    vecs = [x.vec() for x in vecs]
+    M = cls.hstack(*vecs)
+    Q, R = _QRdecomposition_optional(M, normalize=normalize)
 
-        return vec - reduce(lambda a, b: a + b, components)
+    if rankcheck and Q.cols < len(vecs):
+        raise ValueError("GramSchmidt: vector set not linearly independent")
 
-    ret  = []
-    vecs = list(vecs) # make sure we start with a non-zero vector
-
-    while len(vecs) > 0 and vecs[0].is_zero_matrix:
-        if rankcheck is False:
-            del vecs[0]
+    ret = []
+    for i in range(Q.cols):
+        if all_row_vecs:
+            col = cls(Q[:, i].T)
         else:
-            raise ValueError("GramSchmidt: vector set not linearly independent")
-
-    for vec in vecs:
-        perp = perp_to_subspace(vec, ret)
-
-        if not perp.is_zero_matrix:
-            ret.append(cls(perp))
-        elif rankcheck is True:
-            raise ValueError("GramSchmidt: vector set not linearly independent")
-
-    if normalize:
-        ret = [vec / vec.norm() for vec in ret]
-
+            col = cls(Q[:, i])
+        ret.append(col)
     return ret


Test patch:
diff --git a/sympy/matrices/tests/test_decompositions.py b/sympy/matrices/tests/test_decompositions.py
index a39a00c828..383af075c4 100644
--- a/sympy/matrices/tests/test_decompositions.py
+++ b/sympy/matrices/tests/test_decompositions.py
@@ -108,6 +108,12 @@ def test_QR():
     assert R.is_upper
     assert A == Q*R
 
+    A = Matrix([[12, 0, -51], [6, 0, 167], [-4, 0, 24]])
+    Q, R = A.QRdecomposition()
+    assert Q.T * Q == eye(Q.cols)
+    assert R.is_upper
+    assert A == Q*R
+
 def test_QR_non_square():
     # Narrow (cols < rows) matrices
     A = Matrix([[9, 0, 26], [12, 0, -7], [0, 4, 4], [0, -3, -3]])
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 101a16d977..670d72ba53 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -2871,6 +2871,10 @@ def test_gramschmidt_conjugate_dot():
     assert Matrix.orthogonalize(*vecs) == \
         [Matrix([[1], [I]]), Matrix([[1], [-I]])]
 
+    vecs = [Matrix([1, I, 0]), Matrix([I, 0, -I])]
+    assert Matrix.orthogonalize(*vecs) == \
+        [Matrix([[1], [I], [0]]), Matrix([[I/2], [S(1)/2], [-I]])]
+
     mat = Matrix([[1, I], [1, -I]])
     Q, R = mat.QRdecomposition()
     assert Q * Q.H == Matrix.eye(2)