Repository: sympy/sympy
Pull request: Add QR Decomposition Method for DomainMatrix and DDM
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

"Fixes #27250"
#### Brief description of what is fixed or changed
* Added `qr` method to DomainMatrix 
* Implemented Division-based `qr` method for DDM: 
* Added tests to verify:
    * Correct shapes of Q and R
    * Matrix equality A = Q ⋅ R within the specified domain.

#### Other comments
Future work includes adding `qrd` method using fraction free LU decomposition

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* matrices
    * Added qr method for DomainMatrix to compute the QR decomposition.
    * Added a Division-based qr method for DDM to handle domain-aware QR decomposition.
<!-- END RELEASE NOTES -->


Problem statement:
GramSchmidt is really, really slow
I tried applying Gram-Schmidt to construct an orthonormal basis of the following four vectors:

$$
\begin{bmatrix}  -1 + \frac{\sqrt{2}}{2} \\\\  \frac{1}{2} + \frac{\sqrt{2}}{2} \\\\  - \frac{\sqrt{2}}{2} - \frac{1}{2} \\\\  0 \end{bmatrix},
\begin{bmatrix}  -1 - \frac{\sqrt{2}}{2} \\\\  0 \\\\  0 \\\\  0 \end{bmatrix},
\begin{bmatrix}  \frac{\sqrt{2}}{2} \\\\  \frac{1}{2} \\\\  -\frac{1}{2} \\\\  \frac{1}{2} \end{bmatrix},
\begin{bmatrix}  - \frac{\sqrt{2}}{2} \\\\  \frac{\sqrt{2}}{2} \\\\  - \frac{\sqrt{2}}{2} \\\\  -\frac{1}{2} \end{bmatrix}
$$

I've been waiting for _hours_ for a result and don't know why it should take so long.  Here's my code:
```Python
from sympy import *

vecs = [
    Matrix([
        -1 + sqrt(2)/2,
        1/2 + sqrt(2)/2,
        -sqrt(2)/2 - 1/2,
        0,
    ]),
    Matrix([
        -1 - sqrt(2)/2,
        0,
        0,
        0,
    ]),
    Matrix([
        sqrt(2)/2,
        1/2,
        -1/2,
        1/2]),
    Matrix([
        -sqrt(2)/2,
        sqrt(2)/2,
        -sqrt(2)/2,
        -1/2,
    ]),
]
gs = GramSchmidt(vecs, True)
print(gs)
```
Am I doing something wrong, or is this slowness to be expected?  I'm running the SymPy 1.12 package that ships with Ubuntu Noble.

Hints from issue discussion:
It looks like this suffers expression blowup in the QR code. I guess the QR decomposition should be implemented for DomainMatrix.
This is the matrix code for QR adapted to DomainMatrix:
```python
from sympy import Matrix
from sympy.polys.matrices import DomainMatrix
from sympy.polys.matrices.ddm import DDM, DMDomainError


def orthogonalize(A: Matrix):
    Ad: DDM = A.to_DM(extension=True, field=True).to_ddm()

    K = Ad.domain
    if not K.is_Field:
        raise DMDomainError("orthogonalize requires a field")

    rows, cols = Ad.shape
    Q = Ad.copy()
    R = Ad.zeros((rows, cols), K)

    is_zero_col = lambda j: not any(Q[i][j] for i in range(rows))
    dot_cols = lambda i, j: K.sum([Q[k][i] * Q[k][j] for k in range(rows)])

    ranked = []

    for j in range(cols):
        for i in range(j):
            if is_zero_col(i):
                continue

            R[i][j] = dot_cols(i, j) / dot_cols(i, i)
            for k in range(rows):
                Q[k][j] -= R[i][j] * Q[k][i]

        if not is_zero_col(j):
            ranked.append(j)
            R[j][j] = K.one

    Q = Q.extract(range(rows), ranked)
    R = R.extract(ranked, range(cols))

    QM = DomainMatrix.from_rep(Q).to_Matrix()
    RM = DomainMatrix.from_rep(R).to_Matrix()

    return QM, RM
```
Then:
```python
In [2]: M = Matrix.hstack(*vecs).applyfunc(nsimplify)

In [3]: M
Out[3]: 
⎡     √2        √2   √2   -√2 ⎤
⎢-1 + ──   -1 - ──   ──   ────⎥
⎢     2         2    2     2  ⎥
⎢                             ⎥
⎢ 1   √2                   √2 ⎥
⎢ ─ + ──      0     1/2    ── ⎥
⎢ 2   2                    2  ⎥
⎢                             ⎥
⎢  √2   1                 -√2 ⎥
⎢- ── - ─     0     -1/2  ────⎥
⎢  2    2                  2  ⎥
⎢                             ⎥
⎣   0         0     1/2   -1/2⎦

In [4]: edit ortho.py
Editing... done. Executing edited code...

In [5]: %time Q, R = orthogonalize(M)
CPU times: user 24.4 ms, sys: 23 μs, total: 24.4 ms
Wall time: 23.6 ms

In [6]: Q.shape
Out[6]: (4, 3)

In [7]: Q
Out[7]: 
⎡     √2     5   7⋅√2     ⎤
⎢-1 + ──   - ─ - ────   0 ⎥
⎢     2      6    12      ⎥
⎢                         ⎥
⎢ 1   √2     √2   1       ⎥
⎢ ─ + ──   - ── - ──    0 ⎥
⎢ 2   2      12   12      ⎥
⎢                         ⎥
⎢  √2   1   1    √2       ⎥
⎢- ── - ─   ── + ──     0 ⎥
⎢  2    2   12   12       ⎥
⎢                         ⎥
⎣   0          0       1/2⎦

In [8]: expand(Q.T*Q)
Out[8]: 
⎡3     0      0 ⎤
⎢               ⎥
⎢        17     ⎥
⎢0  √2 + ──   0 ⎥
⎢        12     ⎥
⎢               ⎥
⎣0     0     1/4⎦

In [9]: expand(Q*R - M)
Out[9]: 
⎡0  0  0  0⎤
⎢          ⎥
⎢0  0  0  0⎥
⎢          ⎥
⎢0  0  0  0⎥
⎢          ⎥
⎣0  0  0  0⎦
```
There will be some way to adjust the above so that it is fraction-free. I'm not sure what is the best approach to control blowup here:
```python
In [50]: M = randMatrix(3) + randMatrix(3)*x

In [51]: M
Out[51]: 
⎡47⋅x + 56  2⋅x + 94   67⋅x + 88⎤
⎢                               ⎥
⎢88⋅x + 30  68⋅x + 52  11⋅x + 26⎥
⎢                               ⎥
⎣53⋅x + 25  5⋅x + 73   31⋅x + 24⎦

In [52]: %time Q, R = orthogonalize(M)
CPU times: user 65.8 ms, sys: 21 μs, total: 65.8 ms
Wall time: 65.1 ms

In [53]: Q[:,0].T
Out[53]: [47⋅x + 56  88⋅x + 30  53⋅x + 25]

In [54]: Q[:,1].T
Out[54]: 
⎡          3           2                           3          2                               3          2                    ⎤
⎢- 272597⋅x  + 159228⋅x  - 4785⋅x - 46210  309632⋅x  + 38206⋅x  - 212276⋅x - 17098  - 272369⋅x  + 36601⋅x  + 149570⋅x + 124028⎥
⎢────────────────────────────────────────  ───────────────────────────────────────  ──────────────────────────────────────────⎥
⎢              2                                         2                                         2                          ⎥
⎣       12762⋅x  + 13194⋅x + 4661                 12762⋅x  + 13194⋅x + 4661                 12762⋅x  + 13194⋅x + 4661         ⎦

In [55]: Q[:,2].T
Out[55]: 
⎡           5              4               3              2                                     5              4             3              2                                        5              4      ↪
⎢379006068⋅x  + 117583370⋅x  - 1032802148⋅x  + 237364162⋅x  + 265464980⋅x + 31452600  15452523⋅x  - 143100606⋅x  + 25543331⋅x  + 482317556⋅x  - 325550272⋅x - 61420920  - 361756740⋅x  - 105237020⋅x  + 87 ↪
⎢───────────────────────────────────────────────────────────────────────────────────  ────────────────────────────────────────────────────────────────────────────────  ────────────────────────────────── ↪
⎢                     4             3            2                                                        4             3            2                                                       4             ↪
⎣           19147937⋅x  - 26306722⋅x  + 5929529⋅x  - 1204476⋅x + 3821208                        19147937⋅x  - 26306722⋅x  + 5929529⋅x  - 1204476⋅x + 3821208                       19147937⋅x  - 26306722⋅ ↪

↪          3              2                       ⎤
↪ 9870996⋅x  - 349551228⋅x  - 53803192⋅x + 3251280⎥
↪ ────────────────────────────────────────────────⎥
↪  3            2                                 ⎥
↪ x  + 5929529⋅x  - 1204476⋅x + 3821208           ⎦

```
This would be the non-field version:
```python
def orthogonalize(A: Matrix):
    Ad: DDM = A.to_DM(extension=True).to_ddm()

    K = Ad.domain

    rows, cols = Ad.shape
    Q = Ad.copy()

    is_zero_col = lambda j: not any(Q[i][j] for i in range(rows))
    dot_cols = lambda i, j: K.sum([Q[k][i] * Q[k][j] for k in range(rows)])

    ranked = []

    for j in range(cols):
        for i in range(j):
            if is_zero_col(i):
                continue

            _, dij, dii = K.cofactors(dot_cols(i, j), dot_cols(i, i))
            for k in range(rows):
                Q[k][j] = dii * Q[k][j] - dij * Q[k][i]

            g = Q[0][j]
            for k in range(1, rows):
                g = K.gcd(g, Q[k][j])
            for k in range(rows):
                Q[k][j] = K.exquo(Q[k][j], g)

        if not is_zero_col(j):
            ranked.append(j)

    Q = Q.extract(range(rows), ranked)

    QM = DomainMatrix.from_rep(Q).to_Matrix()

    return QM
```
Works better for polynomials but less well for algebraic number fields because gcd there is always 1. There should be some way to generalise this to a fraction-free QR decomposition something like:
```
M = Q*R*D**-1
```
where `D` is a diagonal matrix and `Q`, `R` and `D` are fraction-free.
Fraction-free QR is described here:
https://www.uwo.ca/apmaths/faculty/jeffrey/pdfs/FFLUQR.pdf

> Fraction-free matrix factors: new forms for LU and QR factors
> Wenqin ZHOU (*), David J. JEFFREY
Nice!  Thanks, @oscarbenjamin, for the code.  The various versions do appear to work as advertised.

(I'll have to read up more on domainmatrixes, though.)
Hello @oscarbenjamin , should the `orthogonalize` function replace GramSchmidt, or would implementing QR decomposition for DomainMatrix better address the expression blowup issue you mentioned?
These things should be implemented for DomainMatrix and then the Matrix methods for this should use the DomainMatrix methods.
There should be DomainMatrix methods `qr` and `qrd`. The `qr` method can use the first algorithm I shows and return matrices `Q` and `R` such that `A = Q*R`. The `qrd` method should be more like the second implementation I showed and compute `Q`, `R` and `D` such that `A = Q*D**-1*R` where `D` is a diagonal matrix and `Q` and `R` are in the same domain (integer, polynomial, ...) as the original matrix `A`.

The paper I linked shows how to compute the QRD decomposition using fraction-free LU decomposition which should also be implemented. For now though just adding a `qr` method and ensuring that the matrix code calls it would be good. The implementation I showed can be a method on the `DDM` class that `DomainMatrix` would call.
> There should be DomainMatrix methods `qr` and `qrd`. The `qr` method can use the first algorithm I shows and return matrices `Q` and `R` such that `A = Q*R`. The `qrd` method should be more like the second implementation I showed and compute `Q`, `R` and `D` such that `A = Q*D**-1*R` where `D` is a diagonal matrix and `Q` and `R` are in the same domain (integer, polynomial, ...) as the original matrix `A`.
> 
> The paper I linked shows how to compute the QRD decomposition using fraction-free LU decomposition which should also be implemented. For now though just adding a `qr` method and ensuring that the matrix code calls it would be good. The implementation I showed can be a method on the `DDM` class that `DomainMatrix` would call.

Okay got it, thank you for the clarity

Commit message: Added qr to SDM & DFM

Patch to review:
diff --git a/sympy/polys/matrices/_dfm.py b/sympy/polys/matrices/_dfm.py
index c2f6f16922e0..6a915f4fc61e 100644
--- a/sympy/polys/matrices/_dfm.py
+++ b/sympy/polys/matrices/_dfm.py
@@ -694,6 +694,11 @@ def lu(self):
         L, U, swaps = self.to_ddm().lu()
         return L.to_dfm(), U.to_dfm(), swaps
 
+    def qr(self):
+        """Return the QR decomposition of the matrix."""
+        Q, R = self.to_ddm().qr()
+        return Q.to_dfm(), R.to_dfm()
+
     # XXX: The lu_solve function should be renamed to solve. Whether or not it
     # uses an LU decomposition is an implementation detail. A method called
     # lu_solve would make sense for a situation in which an LU decomposition is
diff --git a/sympy/polys/matrices/ddm.py b/sympy/polys/matrices/ddm.py
index 8f90a3305367..1d027f8a5e6f 100644
--- a/sympy/polys/matrices/ddm.py
+++ b/sympy/polys/matrices/ddm.py
@@ -958,6 +958,34 @@ def lu(a):
 
         return L, U, swaps
 
+    def qr(self):
+        """
+        Fraction-free QR decomposition for DDM.
+        Returns:
+            - Q: Orthogonal matrix as a DDM.
+            - R: Upper triangular matrix as a DDM.
+        """
+        rows, cols = self.shape
+        Q = self.copy()
+        R = self.zeros((rows, cols), self.domain)
+
+        is_zero_col = lambda j: not any(Q[i][j] for i in range(rows))
+        dot_cols = lambda i, j: self.domain.sum(Q[k][i] * Q[k][j] for k in range(rows))
+
+        for j in range(cols):
+            for i in range(j):
+                if is_zero_col(i):
+                    continue
+
+                R[i][j] = dot_cols(i, j) // dot_cols(i, i)
+                for k in range(rows):
+                    Q[k][j] -= R[i][j] * Q[k][i]
+
+            if not is_zero_col(j):
+                R[j][j] = self.domain.one
+
+        return Q, R
+
     def lu_solve(a, b):
         """x where a*x = b"""
         m, n = a.shape
diff --git a/sympy/polys/matrices/domainmatrix.py b/sympy/polys/matrices/domainmatrix.py
index b91bef314d16..ff0bc07248ae 100644
--- a/sympy/polys/matrices/domainmatrix.py
+++ b/sympy/polys/matrices/domainmatrix.py
@@ -3256,6 +3256,23 @@ def lu(self):
         L, U, swaps = self.rep.lu()
         return self.from_rep(L), self.from_rep(U), swaps
 
+    def qr(self):
+        """
+        QR decomposition for DomainMatrix.
+
+        Returns:
+            - Q: Orthogonal DomainMatrix.
+            - R: Upper triangular DomainMatrix.
+        """
+        if isinstance(self.rep, DFM):
+            ddm_q, ddm_r = self.rep.qr()
+        else:
+            ddm_q, ddm_r = self.rep.qr()
+
+        Q = DomainMatrix.from_rep(ddm_q)
+        R = DomainMatrix.from_rep(ddm_r)
+        return Q, R
+
     def lu_solve(self, rhs):
         r"""
         Solver for DomainMatrix x in the A*x = B
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index 0e0685f5d6fb..02aba89e906b 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -1065,6 +1065,19 @@ def lu(A):
         L, U, swaps = A.to_ddm().lu()
         return A.from_ddm(L), A.from_ddm(U), swaps
 
+    def qr(self):
+        """
+        QR decomposition for SDM (Sparse Domain Matrix).
+
+        Returns:
+            - Q: Orthogonal matrix as a SDM.
+            - R: Upper triangular matrix as a SDM.
+        """
+        ddm_q, ddm_r = self.to_ddm().qr()
+        Q = ddm_q.to_sdm()
+        R = ddm_r.to_sdm()
+        return Q, R
+
     def lu_solve(A, b):
         """
 
diff --git a/sympy/polys/matrices/tests/test_ddm.py b/sympy/polys/matrices/tests/test_ddm.py
index 44c862461e85..5bb1ba920467 100644
--- a/sympy/polys/matrices/tests/test_ddm.py
+++ b/sympy/polys/matrices/tests/test_ddm.py
@@ -7,6 +7,7 @@
 from sympy.polys.matrices.exceptions import (
     DMShapeError, DMNonInvertibleMatrixError, DMDomainError,
     DMBadInputError)
+from sympy.polys.matrices import DomainMatrix
 
 
 def test_DDM_init():
@@ -556,3 +557,15 @@ def test_DDM_is_lower():
     ], (4, 3), QQ).transpose()
     assert A.is_lower() is True
     assert B.is_lower() is False
+
+
+def test_ddm_qr():
+    # Create a sample matrix in DomainMatrix
+    A = DomainMatrix([[ZZ(3), ZZ(1)], [ZZ(4), ZZ(3)]], (2, 2), ZZ)
+    # Perform QR decomposition
+    Q, R = A.qr()
+    # Check that the shapes are correct
+    assert Q.shape == (2, 2)
+    assert R.shape == (2, 2)
+    # Check that A = Q * R (within the domain of integers)
+    assert Q * R == A
diff --git a/sympy/polys/matrices/tests/test_domainmatrix.py b/sympy/polys/matrices/tests/test_domainmatrix.py
index 2b59d76a9d46..cd305f492eb8 100644
--- a/sympy/polys/matrices/tests/test_domainmatrix.py
+++ b/sympy/polys/matrices/tests/test_domainmatrix.py
@@ -1360,3 +1360,11 @@ def test_DomainMatrix_pickling():
     assert pickle.loads(pickle.dumps(dM)) == dM
     dM = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
     assert pickle.loads(pickle.dumps(dM)) == dM
+
+
+def test_qr():
+    A = DomainMatrix([[ZZ(3), ZZ(1)], [ZZ(4), ZZ(3)]], (2, 2), ZZ)
+    Q, R = A.qr()
+    assert Q.shape == (2, 2)
+    assert R.shape == (2, 2)
+    assert Q * R == A


Code review comments:
diff_hunk: @@ -3256,6 +3256,23 @@ def lu(self):
         L, U, swaps = self.rep.lu()
         return self.from_rep(L), self.from_rep(U), swaps
 
+    def qr(self):
+        """
+        QR decomposition for DomainMatrix.
+
+        Returns:
+            - Q: Orthogonal DomainMatrix.
+            - R: Upper triangular DomainMatrix.
+        """
+        if isinstance(self.rep, DFM):
+            ddm_q, ddm_r = self.rep.qr()
+        else:
+            ddm_q, ddm_r = self.rep.qr()
path: sympy/polys/matrices/domainmatrix.py
text: @user1:
If the code is the same in both branches then the if check is not needed.

@author:
Okay, I will update that.

diff_hunk: @@ -958,6 +958,34 @@ def lu(a):
 
         return L, U, swaps
 
+    def qr(self):
+        """
+        Fraction-free QR decomposition for DDM.
+        Returns:
+            - Q: Orthogonal matrix as a DDM.
+            - R: Upper triangular matrix as a DDM.
+        """
+        rows, cols = self.shape
+        Q = self.copy()
+        R = self.zeros((rows, cols), self.domain)
+
+        is_zero_col = lambda j: not any(Q[i][j] for i in range(rows))
+        dot_cols = lambda i, j: self.domain.sum(Q[k][i] * Q[k][j] for k in range(rows))
+
+        for j in range(cols):
+            for i in range(j):
+                if is_zero_col(i):
+                    continue
+
+                R[i][j] = dot_cols(i, j) // dot_cols(i, i)
path: sympy/polys/matrices/ddm.py
text: @user1:
Using `//` here is incorrect. This is not how fraction-free algorithms work: exact division (`exquo`) can be used but floor division should not.

@author:
Okay, I used `/` at first but after installing `python-flint`, `/` started giving  me  this` "error:  flint.utils.flint_exceptions.DomainError: fmpz division is not exact "`


@author:
but i will work around it, thank you for the nuggets. 

@user1:
If the domain is a field then `/` can be used.

@author:
what are the domains the `qr` method should support?  The current implementation only supports integers (ZZ).

@user1:
The QR decomposition should be computable for any field. To learn more about the domains see:
https://docs.sympy.org/latest/modules/polys/domainsintro.html