'Return Q, R where A = Q*R, Q is orthogonal and R is upper triangular.
Examples
This is the example from wikipedia:
>>> from sympy import Matrix
>>> A = Matrix([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])
>>> Q, R = A.QRdecomposition()
>>> Q
Matrix([
[ 6/7, -69/175, -58/175],
[ 3/7, 158/175,   6/175],
[-2/7,    6/35,  -33/35]])
>>> R
Matrix([
[14,  21, -14],
[ 0, 175, -70],
[ 0,   0,  35]])
>>> A == Q*R
True
QR factorization of an identity matrix:
>>> A = Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
>>> Q, R = A.QRdecomposition()
>>> Q
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> R
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
See Also
cholesky
LDLdecomposition
LUdecomposition
QRsolve'

def QRdecomposition(self):

    cls = self.__class__
    mat = self.as_mutable()
    if (not (mat.rows >= mat.cols)):
        raise MatrixError('The    number    of    rows    must    be    greater    than    columns')
    n = mat.rows
    m = mat.cols
    rank = n
    row_reduced = mat.rref()[0]
    for i in range(row_reduced.rows):
        if (row_reduced.row(i).norm() == 0):
            rank -= 1
    if (not (rank == mat.cols)):
        raise MatrixError('The    rank    of    the    matrix    must    match    the    columns')
    (Q, R) = (mat.zeros(n, m), mat.zeros(m))
    for j in range(m):
        tmp = mat[:, j]
        for i in range(j):
            tmp -= (Q[:, i] * mat[:, j].dot(Q[:, i]))
            tmp.expand()
        R[(j, j)] = tmp.norm()
        Q[:, j] = (tmp / R[(j, j)])
        if (Q[:, j].norm() != 1):
            raise NotImplementedError(('Could    not    normalize    the    vector    %d.' % j))
        for i in range(j):
            R[(i, j)] = Q[:, i].dot(mat[:, j])
    return (cls(Q), cls(R))
