MAT 232 Chapter 6

76 days ago by kcrisman

Chapter 6 - Orthogonality

Section 6.1

There are many ways for things to be orthogonal.

vector([1,2,3]).dot_product(vector([-3,0,1])) 
       
0
0
vector([1,2,3]) * vector([-3,0,1]) 
       
0
0

For instance, we have some polynomials that are, if you consider the inner product to be $\int_{-1}^1 f(x)g(x)dx$.

for i in range(4): show(expand(legendre_P(i,x))) 
       



for i in range(4): for j in range(i+1,4): f(x) = legendre_P(i,x) g(x) = legendre_P(j,x) html("$$\int_{-1}^1 \left(%s\\right) \left(%s\\right)\; dx = %s$$"%(latex(f(x)),latex(g(x)),integrate(legendre_P(i,x)*legendre_P(j,x),(x,-1,1)))) 
       




















Or you can try some complex vectors.

vector([1+i,2])*vector([1-i,4]) 
       
10
10
vector([1,i]) * vector([1,i]) 
       
0
0

The functions $\sin(kx)$ and $\cos(\ell x)$ are also all orthogonal on $[0,2\pi]$.

@interact def _(i=([1..5]),j=([1..5])): html("The integral $$\int_0^{2\pi} \sin(%s x)\sin(%s x)\; dx$$"%(i,j)) html("is $%s$"%latex(integrate(sin(i*x)*sin(j*x),(x,0,2*pi)))) html("The integral $$\int_0^{2\pi} \cos(%s x)\sin(%s x)\; dx$$"%(i,j)) html("is $%s$"%latex(integrate(cos(i*x)*sin(j*x),(x,0,2*pi)))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Anyway, now we have a way to talk about what we claimed before - that the row space and null space are "perpendicular"!

M = matrix([[1,2,3],[1,2,3],[1,2,3]]); show(M) 
       
V = M.row_space() W = M.right_kernel() 
       
V;W 
       
Free module of degree 3 and rank 1 over Integer Ring
Echelon basis matrix:
[1 2 3]
Free module of degree 3 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1  1 -1]
[ 0  3 -2]
Free module of degree 3 and rank 1 over Integer Ring
Echelon basis matrix:
[1 2 3]
Free module of degree 3 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1  1 -1]
[ 0  3 -2]

We just multiply every basis vector of the row space with every basis vector of the null space.

for v in V.basis(): for w in W.basis(): html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w))) 
       




@interact def _(M = input_grid(4, 5, default=[[-1,-2,1,-1,0],[3,6,6,0,3],[3,6,4,1,3],[-2,-4,-5,0,-3]])): M = matrix(M) V = M.right_kernel() W = M.row_space() for v in V.basis(): for w in W.basis(): pretty_print(html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w)))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Incidentally, the same is true of the transpose (the column space and kernel of the transpose):

@interact def _(M = input_grid(4, 5, default=[[-1,-2,1,-1,0],[3,6,6,0,3],[3,6,4,1,3],[-2,-4,-5,0,-3]])): M = matrix(M) V = M.left_kernel() W = M.column_space() html("Now we're doing it with the column space and <em>left</em> kernel") for v in V.basis(): for w in W.basis(): html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Warning!  There is a lot we are not saying about these complements - see Sage ticket #7522 for some discussion about this between those interested in undergraduate pedagogical situations and those warning of keeping things completely mathematically correct!

Section 6.2

Let's do some orthogonal projection.

@interact def _(V = matrix([[1,3]]),W = matrix([[-2,3]])): V = vector(V.list()) W = vector(W.list()) ProjVW = (V*W)/(W*W)*W show(plot(ProjVW,color='green',zorder=5)+plot(V,color='green')+plot(W)+line([V.list(),ProjVW.list()],color='black',linestyle='--'),aspect_ratio=1) 
       

Click to the left again to hide and once more to show the dynamic interactive window

 
       

Section 6.3

Well, what about in more than two dimensions?

@interact def _(V = matrix([[1,3,2]]),W = matrix([[-2,3,2]])): V = vector(V.list()) W = vector(W.list()) ProjVW = (V*W)/(W*W)*W show(plot(ProjVW,color='green',zorder=5)+plot(V,color='green')+plot(W)+line([V.list(),ProjVW.list()],color='black',linestyle='--'),aspect_ratio=1) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Note that the default setting here doesn't have the $W1$ and $W2$ being orthogonal!  You can't use the formula without doing that first.

var('u,v') @interact def _(V = matrix([[1,3,2]]),W1 = matrix([[-2,3,2]]),W2 = matrix([[-2,1,1]])): V = vector(V.list()) W1 = vector(W1.list()) W2 = vector(W2.list()) if W1*W2 != 0: print "you don't have orthogonal vectors" else: ProjVW = (V*W1)/(W1*W1)*W1+(V*W2)/(W2*W2)*W2 show(plot(ProjVW,color='green',zorder=5)+plot(V,color='blue',thickness=2)+parametric_plot3d(u*W1+v*W2, (u,-3,3), (v,-3,3),color='green',opacity=.7)+line([V.list(),ProjVW.list()],color='black',linestyle='--',thickness=5),aspect_ratio=1) 
       

Click to the left again to hide and once more to show the dynamic interactive window

Of course, we can do all this for the polynomial or other inner products...

w1 = x w2 = 3/2*x^2-1/2 v = 1+x+x^2 
       
F = integrate(w1*v,x,-1,1)/integrate(w1*w1,x,-1,1)*w1+integrate(w2*v,x,-1,1)/integrate(w2*w2,x,-1,1)*w2 F 
       
x^2 + x - 1/3
x^2 + x - 1/3

Here, $F$ is the projection and $v-F$ is the residual.

integrate((x^2+x-1/3)*(4/3),x,-1,1) 
       
0
0
integrate((v-F)*F,x,-1,1) 
       
0
0

Inner product is zero.  Yeah.  And the distance from $v$ to this subspace is...

sqrt(integrate((v-F)^2,x,-1,1)) 
       
4/3*sqrt(2)
4/3*sqrt(2)

Huh.

Section 6.4

So, how does one exactly get this orthogonal basis?  The answer is the Gram-Schmidt process!  

#A = matrix([[1,1,0,1,0],[1,0,1,1,0]]) #A = matrix([[1,1,0,1,0],[1,0,1,1,0],[0,1,1,1,0]]) A = matrix([[1,1,0,1,0],[1,0,1,1,0],[0,1,1,1,0],[5,0,0,0,1]]) show(A) 
       

Sage returns the transpose of the column vectors the book might lead us to expect.

show(A.gram_schmidt()[0].transpose()) 
       

And Sage actually returns a matrix of orthogonal rows along with the lower triangular matrix of the transformations needed to achieve it.

show(A.gram_schmidt()) 
       

Let's see what this looks like for three general vectors.

var('u,v') @interact def _(U1 = matrix([[1,3,2]]),U2 = matrix([[-2,3,2]]),U3 = matrix([[-2,1,1]])): U1 = vector(U1.list()) U2 = vector(U2.list()) U3 = vector(U3.list()) M = matrix([U1,U2,U3]) if M.det()==0: print "you don't have a basis!" else: M = M.gram_schmidt()[0] G = Graphics() G += plot(U1,thickness=2)+plot(U2,thickness=2)+plot(U3,thickness=2) G += parametric_plot3d(u*U1+v*U2, (u,-1,1), (v,-1,1),color='green',opacity=.6) V1,V2,V3 = M.rows() G += plot(V1,thickness=2,color='red')+plot(V2,thickness=2,color='red')+plot(V3,thickness=2,color='red') G += line([V2.list(),U2.list()],color='black',linestyle='--',thickness=2) G += line([V3.list(),U3.list()],color='black',linestyle='--',thickness=2) show(G,aspect_ratio=1) 
       

Click to the left again to hide and once more to show the dynamic interactive window

We can use a similar method to get a set of orthogonal polynomials - very useful in applications.

legendre_P? 
       

File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.py

Type: <type ‘function’>

Definition: legendre_P(n, x)

Docstring:

Returns the Legendre polynomial of the first kind for integers \(n > -1\).

REFERENCE:

  • AS 22.5.35 page 779.

EXAMPLES:

sage: P.<t> = QQ[]
sage: legendre_P(2,t)
3/2*t^2 - 1/2
sage: legendre_P(3, 1.1)
1.67750000000000
sage: legendre_P(3, 1 + I)
7/2*I - 13/2
sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
[-179  242]
[-484  547]
sage: legendre_P(3, GF(11)(5))
8

File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.py

Type: <type ‘function’>

Definition: legendre_P(n, x)

Docstring:

Returns the Legendre polynomial of the first kind for integers \(n > -1\).

REFERENCE:

  • AS 22.5.35 page 779.

EXAMPLES:

sage: P.<t> = QQ[]
sage: legendre_P(2,t)
3/2*t^2 - 1/2
sage: legendre_P(3, 1.1)
1.67750000000000
sage: legendre_P(3, 1 + I)
7/2*I - 13/2
sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
[-179  242]
[-484  547]
sage: legendre_P(3, GF(11)(5))
8

Sage also has the QR decomposition for numerical matrices.

show(A.change_ring(RDF)) show(A.change_ring(RDF).QR()) 
       

Note that it returns something similar to Matlab, not quite what the book has; see Exercise 21 in your text.

A.QR? 
       

File: /usr/local/sage-6.9/src/sage/matrix/matrix2.pyx

Type: <type ‘builtin_function_or_method’>

Definition: A.QR(full=True)

Docstring:

Returns a factorization of self as a unitary matrix and an upper-triangular matrix.

INPUT:

  • full - default: True - if True then the returned matrices have dimensions as described below. If False the R matrix has no zero rows and the columns of Q are a basis for the column space of self.

OUTPUT:

If self is an m\times n matrix and full=True then this method returns a pair of matrices: Q is an m\times m unitary matrix (meaning its inverse is its conjugate-transpose) and R is an m\times n upper-triangular matrix with non-negative entries on the diagonal. For a matrix of full rank this factorization is unique (due to the restriction to positive entries on the diagonal).

If full=False then Q has m rows and the columns form an orthonormal basis for the column space of self. So, in particular, the conjugate-transpose of Q times Q will be an identity matrix. The matrix R will still be upper-triangular but will also have full rank, in particular it will lack the zero rows present in a full factorization of a rank-deficient matrix.

The results obtained when full=True are cached, hence Q and R are immutable matrices in this case.

Note

This is an exact computation, so limited to exact rings. Also the base ring needs to have a fraction field implemented in Sage and this field must contain square roots. One example is the field of algebraic numbers, QQbar, as used in the examples below. If you need numerical results, convert the base ring to the field of complex double numbers, CDF, which will use a faster routine that is careful about numerical subtleties.

ALGORITHM:

“Modified Gram-Schmidt,” Algorithm 8.1 of [TREFETHEN-BAU].

EXAMPLES:

For a nonsingular matrix, the QR decomposition is unique.

sage: A = matrix(QQbar, [[-2, 0, -4, -1, -1],
...                      [-2, 1, -6, -3, -1],
...                      [1, 1, 7, 4, 5],
...                      [3, 0, 8, 3, 3],
...                      [-1, 1, -6, -6, 5]])
sage: Q, R = A.QR()
sage: Q
[ -0.4588314677411235?  -0.1260506983326509?   0.3812120831224489?   -0.394573711338418?      -0.687440062597?]
[ -0.4588314677411235?   0.4726901187474409? -0.05198346588033394?    0.717294125164660?      -0.220962877263?]
[  0.2294157338705618?   0.6617661662464172?   0.6619227988762521?   -0.180872093737548?      0.1964114464561?]
[  0.6882472016116853?   0.1890760474989764?  -0.2044682991293135?    0.096630296654307?      -0.662888631790?]
[ -0.2294157338705618?   0.5357154679137663?   -0.609939332995919?   -0.536422031427112?      0.0245514308070?]
sage: R
[  4.358898943540674? -0.4588314677411235?   13.07669683062202?   6.194224814505168?   2.982404540317303?]
[                   0   1.670171752907625?  0.5987408170800917?  -1.292019657909672?   6.207996892883057?]
[                   0                    0   5.444401659866974?   5.468660610611130?  -0.682716185228386?]
[                   0                    0                    0   1.027626039419836?   -3.61930014968662?]
[                   0                    0                    0                    0    0.02455143080702?]
sage: Q.conjugate_transpose()*Q
[1.000000000000000?            0.?e-18            0.?e-17            0.?e-15            0.?e-12]
[           0.?e-18 1.000000000000000?            0.?e-16            0.?e-15            0.?e-12]
[           0.?e-17            0.?e-16 1.000000000000000?            0.?e-15            0.?e-12]
[           0.?e-15            0.?e-15            0.?e-15 1.000000000000000?            0.?e-12]
[           0.?e-12            0.?e-12            0.?e-12            0.?e-12    1.000000000000?]
sage: Q*R == A
True

An example with complex numbers in QQbar, the field of algebraic numbers.

sage: A = matrix(QQbar, [[-8, 4*I + 1, -I + 2, 2*I + 1],
...                      [1, -2*I - 1, -I + 3, -I + 1],
...                      [I + 7, 2*I + 1, -2*I + 7, -I + 1],
...                      [I + 2, 0, I + 12, -1]])
sage: Q, R = A.QR()
sage: Q
[                          -0.7302967433402215?    0.2070566455055649? + 0.5383472783144687?*I    0.2463049809998642? - 0.0764456358723292?*I    0.2381617683194332? - 0.1036596032779695?*I]
[                           0.0912870929175277?   -0.2070566455055649? - 0.3778783780476559?*I    0.3786559533863032? - 0.1952221495524667?*I      0.701244450214469? - 0.364371165098660?*I]
[   0.6390096504226938? + 0.0912870929175277?*I    0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I    0.3140171085506763? - 0.0825191718705412?*I]
[   0.1825741858350554? + 0.0912870929175277?*I  -0.03623491296347385? + 0.0724698259269477?*I   0.8632284069415110? + 0.06322839976356195?*I   -0.4499694867611521? - 0.0116119181208918?*I]
sage: R
[                          10.95445115010333?               0.?e-18 - 1.917028951268082?*I    5.385938482134133? - 2.190890230020665?*I  -0.2738612787525831? - 2.190890230020665?*I]
[                                           0               4.829596256417300? + 0.?e-17*I   -0.869637911123373? - 5.864879483945125?*I   0.993871898426712? - 0.3054085521207082?*I]
[                                           0                                            0               12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I]
[                                           0                                            0                                            0               1.942963944258992? + 0.?e-16*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-19*I            0.?e-18 + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-17 + 0.?e-17*I            0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-16*I            0.?e-16 + 0.?e-16*I]
[           0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-15*I]
sage: Q*R - A
[            0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[            0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-15*I]
[0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-16*I]

A rank-deficient rectangular matrix, with both values of the full keyword.

sage: A = matrix(QQbar, [[2, -3, 3],
...                      [-1, 1, -1],
...                      [-1, 3, -3],
...                      [-5, 1, -1]])
sage: Q, R = A.QR()
sage: Q
[  0.3592106040535498?  -0.5693261797050169?   0.7239227659930268?   0.1509015305256380?]
[ -0.1796053020267749?   0.1445907757980996?                     0   0.9730546968377341?]
[ -0.1796053020267749?   0.7048800320157352?    0.672213996993525?  -0.1378927778941174?]
[ -0.8980265101338745?  -0.3976246334447737?   0.1551263069985058? -0.10667177157846818?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
[                  0                   0                   0]
[                  0                   0                   0]
sage: Q.conjugate_transpose()*Q
[                 1            0.?e-18            0.?e-18            0.?e-18]
[           0.?e-18                  1            0.?e-18            0.?e-18]
[           0.?e-18            0.?e-18 1.000000000000000?            0.?e-18]
[           0.?e-18            0.?e-18            0.?e-18 1.000000000000000?]

sage: Q, R = A.QR(full=False)
sage: Q
[ 0.3592106040535498? -0.5693261797050169?]
[-0.1796053020267749?  0.1445907757980996?]
[-0.1796053020267749?  0.7048800320157352?]
[-0.8980265101338745? -0.3976246334447737?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
sage: Q.conjugate_transpose()*Q
[      1 0.?e-18]
[0.?e-18       1]

Another rank-deficient rectangular matrix, with complex entries, as a reduced decomposition.

sage: A = matrix(QQbar, [[-3*I - 3, I - 3, -12*I + 1, -2],
...                      [-I - 1, -2, 5*I - 1, -I - 2],
...                      [-4*I - 4, I - 5, -7*I, -I - 4]])
sage: Q, R = A.QR(full=False)
sage: Q
[ -0.4160251471689219? - 0.4160251471689219?*I   0.5370861555295747? + 0.1790287185098583?*I]
[ -0.1386750490563073? - 0.1386750490563073?*I  -0.7519206177414046? - 0.2506402059138015?*I]
[ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I]
sage: R
[                        7.211102550927979?  3.328201177351375? - 5.269651864139676?*I   7.904477796209515? + 8.45917799243475?*I  4.021576422632911? - 2.634825932069838?*I]
[                                         0                         1.074172311059150?  -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-18*I            0.?e-18 + 0.?e-18*I]
[           0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I]
sage: Q*R-A
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]
[0.?e-18 + 0.?e-18*I 0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]

Results of full decompositions are cached and thus returned immutable.

sage: A = random_matrix(QQbar, 2, 2)
sage: Q, R = A.QR()
sage: Q.is_mutable()
False
sage: R.is_mutable()
False

Trivial cases return trivial results of the correct size, and we check Q itself in one case.

sage: A = zero_matrix(QQbar, 0, 10)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(0, 0)
sage: R.nrows(), R.ncols()
(0, 10)
sage: A = zero_matrix(QQbar, 3, 0)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(3, 3)
sage: R.nrows(), R.ncols()
(3, 0)
sage: Q
[1 0 0]
[0 1 0]
[0 0 1]

TESTS:

Inexact rings are caught and CDF suggested.

sage: A = matrix(RealField(100), 2, range(4))
sage: A.QR()
Traceback (click to the left of this block for traceback)
...
                                
                            

File: /usr/local/sage-6.9/src/sage/matrix/matrix2.pyx

Type: <type ‘builtin_function_or_method’>

Definition: A.QR(full=True)

Docstring:

Returns a factorization of self as a unitary matrix and an upper-triangular matrix.

INPUT:

  • full - default: True - if True then the returned matrices have dimensions as described below. If False the R matrix has no zero rows and the columns of Q are a basis for the column space of self.

OUTPUT:

If self is an m\times n matrix and full=True then this method returns a pair of matrices: Q is an m\times m unitary matrix (meaning its inverse is its conjugate-transpose) and R is an m\times n upper-triangular matrix with non-negative entries on the diagonal. For a matrix of full rank this factorization is unique (due to the restriction to positive entries on the diagonal).

If full=False then Q has m rows and the columns form an orthonormal basis for the column space of self. So, in particular, the conjugate-transpose of Q times Q will be an identity matrix. The matrix R will still be upper-triangular but will also have full rank, in particular it will lack the zero rows present in a full factorization of a rank-deficient matrix.

The results obtained when full=True are cached, hence Q and R are immutable matrices in this case.

Note

This is an exact computation, so limited to exact rings. Also the base ring needs to have a fraction field implemented in Sage and this field must contain square roots. One example is the field of algebraic numbers, QQbar, as used in the examples below. If you need numerical results, convert the base ring to the field of complex double numbers, CDF, which will use a faster routine that is careful about numerical subtleties.

ALGORITHM:

“Modified Gram-Schmidt,” Algorithm 8.1 of [TREFETHEN-BAU].

EXAMPLES:

For a nonsingular matrix, the QR decomposition is unique.

sage: A = matrix(QQbar, [[-2, 0, -4, -1, -1],
...                      [-2, 1, -6, -3, -1],
...                      [1, 1, 7, 4, 5],
...                      [3, 0, 8, 3, 3],
...                      [-1, 1, -6, -6, 5]])
sage: Q, R = A.QR()
sage: Q
[ -0.4588314677411235?  -0.1260506983326509?   0.3812120831224489?   -0.394573711338418?      -0.687440062597?]
[ -0.4588314677411235?   0.4726901187474409? -0.05198346588033394?    0.717294125164660?      -0.220962877263?]
[  0.2294157338705618?   0.6617661662464172?   0.6619227988762521?   -0.180872093737548?      0.1964114464561?]
[  0.6882472016116853?   0.1890760474989764?  -0.2044682991293135?    0.096630296654307?      -0.662888631790?]
[ -0.2294157338705618?   0.5357154679137663?   -0.609939332995919?   -0.536422031427112?      0.0245514308070?]
sage: R
[  4.358898943540674? -0.4588314677411235?   13.07669683062202?   6.194224814505168?   2.982404540317303?]
[                   0   1.670171752907625?  0.5987408170800917?  -1.292019657909672?   6.207996892883057?]
[                   0                    0   5.444401659866974?   5.468660610611130?  -0.682716185228386?]
[                   0                    0                    0   1.027626039419836?   -3.61930014968662?]
[                   0                    0                    0                    0    0.02455143080702?]
sage: Q.conjugate_transpose()*Q
[1.000000000000000?            0.?e-18            0.?e-17            0.?e-15            0.?e-12]
[           0.?e-18 1.000000000000000?            0.?e-16            0.?e-15            0.?e-12]
[           0.?e-17            0.?e-16 1.000000000000000?            0.?e-15            0.?e-12]
[           0.?e-15            0.?e-15            0.?e-15 1.000000000000000?            0.?e-12]
[           0.?e-12            0.?e-12            0.?e-12            0.?e-12    1.000000000000?]
sage: Q*R == A
True

An example with complex numbers in QQbar, the field of algebraic numbers.

sage: A = matrix(QQbar, [[-8, 4*I + 1, -I + 2, 2*I + 1],
...                      [1, -2*I - 1, -I + 3, -I + 1],
...                      [I + 7, 2*I + 1, -2*I + 7, -I + 1],
...                      [I + 2, 0, I + 12, -1]])
sage: Q, R = A.QR()
sage: Q
[                          -0.7302967433402215?    0.2070566455055649? + 0.5383472783144687?*I    0.2463049809998642? - 0.0764456358723292?*I    0.2381617683194332? - 0.1036596032779695?*I]
[                           0.0912870929175277?   -0.2070566455055649? - 0.3778783780476559?*I    0.3786559533863032? - 0.1952221495524667?*I      0.701244450214469? - 0.364371165098660?*I]
[   0.6390096504226938? + 0.0912870929175277?*I    0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I    0.3140171085506763? - 0.0825191718705412?*I]
[   0.1825741858350554? + 0.0912870929175277?*I  -0.03623491296347385? + 0.0724698259269477?*I   0.8632284069415110? + 0.06322839976356195?*I   -0.4499694867611521? - 0.0116119181208918?*I]
sage: R
[                          10.95445115010333?               0.?e-18 - 1.917028951268082?*I    5.385938482134133? - 2.190890230020665?*I  -0.2738612787525831? - 2.190890230020665?*I]
[                                           0               4.829596256417300? + 0.?e-17*I   -0.869637911123373? - 5.864879483945125?*I   0.993871898426712? - 0.3054085521207082?*I]
[                                           0                                            0               12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I]
[                                           0                                            0                                            0               1.942963944258992? + 0.?e-16*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-19*I            0.?e-18 + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I            0.?e-17 + 0.?e-17*I            0.?e-16 + 0.?e-16*I]
[           0.?e-17 + 0.?e-17*I            0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-16*I            0.?e-16 + 0.?e-16*I]
[           0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I            0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-15*I]
sage: Q*R - A
[            0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[            0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-15*I]
[0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-16*I]

A rank-deficient rectangular matrix, with both values of the full keyword.

sage: A = matrix(QQbar, [[2, -3, 3],
...                      [-1, 1, -1],
...                      [-1, 3, -3],
...                      [-5, 1, -1]])
sage: Q, R = A.QR()
sage: Q
[  0.3592106040535498?  -0.5693261797050169?   0.7239227659930268?   0.1509015305256380?]
[ -0.1796053020267749?   0.1445907757980996?                     0   0.9730546968377341?]
[ -0.1796053020267749?   0.7048800320157352?    0.672213996993525?  -0.1378927778941174?]
[ -0.8980265101338745?  -0.3976246334447737?   0.1551263069985058? -0.10667177157846818?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
[                  0                   0                   0]
[                  0                   0                   0]
sage: Q.conjugate_transpose()*Q
[                 1            0.?e-18            0.?e-18            0.?e-18]
[           0.?e-18                  1            0.?e-18            0.?e-18]
[           0.?e-18            0.?e-18 1.000000000000000?            0.?e-18]
[           0.?e-18            0.?e-18            0.?e-18 1.000000000000000?]

sage: Q, R = A.QR(full=False)
sage: Q
[ 0.3592106040535498? -0.5693261797050169?]
[-0.1796053020267749?  0.1445907757980996?]
[-0.1796053020267749?  0.7048800320157352?]
[-0.8980265101338745? -0.3976246334447737?]
sage: R
[ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
[                  0  3.569584777515583? -3.569584777515583?]
sage: Q.conjugate_transpose()*Q
[      1 0.?e-18]
[0.?e-18       1]

Another rank-deficient rectangular matrix, with complex entries, as a reduced decomposition.

sage: A = matrix(QQbar, [[-3*I - 3, I - 3, -12*I + 1, -2],
...                      [-I - 1, -2, 5*I - 1, -I - 2],
...                      [-4*I - 4, I - 5, -7*I, -I - 4]])
sage: Q, R = A.QR(full=False)
sage: Q
[ -0.4160251471689219? - 0.4160251471689219?*I   0.5370861555295747? + 0.1790287185098583?*I]
[ -0.1386750490563073? - 0.1386750490563073?*I  -0.7519206177414046? - 0.2506402059138015?*I]
[ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I]
sage: R
[                        7.211102550927979?  3.328201177351375? - 5.269651864139676?*I   7.904477796209515? + 8.45917799243475?*I  4.021576422632911? - 2.634825932069838?*I]
[                                         0                         1.074172311059150?  -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-18*I            0.?e-18 + 0.?e-18*I]
[           0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I]
sage: Q*R-A
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]
[0.?e-18 + 0.?e-18*I 0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-18 + 0.?e-18*I]

Results of full decompositions are cached and thus returned immutable.

sage: A = random_matrix(QQbar, 2, 2)
sage: Q, R = A.QR()
sage: Q.is_mutable()
False
sage: R.is_mutable()
False

Trivial cases return trivial results of the correct size, and we check Q itself in one case.

sage: A = zero_matrix(QQbar, 0, 10)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(0, 0)
sage: R.nrows(), R.ncols()
(0, 10)
sage: A = zero_matrix(QQbar, 3, 0)
sage: Q, R = A.QR()
sage: Q.nrows(), Q.ncols()
(3, 3)
sage: R.nrows(), R.ncols()
(3, 0)
sage: Q
[1 0 0]
[0 1 0]
[0 0 1]

TESTS:

Inexact rings are caught and CDF suggested.

sage: A = matrix(RealField(100), 2, range(4))
sage: A.QR()
Traceback (most recent call last):
...
NotImplementedError: QR decomposition is implemented over exact rings, try CDF for numerical results, not Real Field with 100 bits of precision

Without a fraction field, we cannot hope to run the algorithm.

sage: A = matrix(Integers(6), 2, range(4))
sage: A.QR()
Traceback (most recent call last):
...
ValueError: QR decomposition needs a fraction field of Ring of integers modulo 6

The biggest obstacle is making unit vectors, thus requiring square roots, though some small cases pass through.

sage: A = matrix(ZZ, 3, range(9))
sage: A.QR()
Traceback (most recent call last):
...
TypeError: QR decomposition unable to compute square roots in Rational Field

sage: A = matrix(ZZ, 2, range(4))
sage: Q, R = A.QR()
sage: Q
[0 1]
[1 0]
sage: R
[2 3]
[0 1]

REFERENCES:

[TREFETHEN-BAU]Trefethen, Lloyd N., Bau, David, III “Numerical Linear Algebra” SIAM, Philadelphia, 1997.

AUTHOR:

  • Rob Beezer (2011-02-17)

But not symbolic or integer ones.

A.QR() 
       
Traceback (click to the left of this block for traceback)
...
TypeError: QR decomposition unable to compute square roots in Rational
Field
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_19.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QS5RUigp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmplmJ82W/___code___.py", line 2, in <module>
    exec compile(u'A.QR()
  File "", line 1, in <module>
    
  File "sage/matrix/matrix2.pyx", line 8978, in sage.matrix.matrix2.Matrix.QR (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix2.c:65263)
TypeError: QR decomposition unable to compute square roots in Rational Field

Sections 6.5/6.6

Let's fit some linear models, implicitly using least squares.

var('a,b') model(x) = a+b*x Data = [(1,0),(2,1),(4,2),(5,3)] find_fit(Data,model) 
       
[a == -0.6000000000034891, b == 0.6999999999993456]
[a == -0.6000000000034891, b == 0.6999999999993456]
M = find_fit(Data,model) points(Data)+plot(model.subs(a=M[0].rhs(),b=M[1].rhs()),(x,0,5),color='black') 
       

Here's some real-life data about men's and women's record times in the 100-meter free in swimming.

var('a,b,c') model(gender,year) = a+b*year+c*gender Data = [(1,1905,65.8),(1,1908,65.6),(1,1910,62.8),(1,1912,61.6),(1,1918,61.4),(1,1920,60.4),(1,1922,58.6),(1,1924,57.4),(1,1934,56.8),(1,1935,56.6),(1,1936,56.4),(1,1944,55.9),(1,1947,55.8),(1,1948,55.4),(1,1955,54.8),(1,1957,54.6),(1,1961,53.6),(1,1964,52.9),(1,1967,52.6),(1,1968,52.2),(1,1970,51.9),(1,1972,51.22),(1,1975,50.59),(1,1976,49.44),(1,1981,49.36),(1,1985,49.24),(1,1986,48.74),(1,1988,48.42),(1,1994,48.21),(1,2000,48.18),(1,2000,47.84),(0,1908,95),(0,1910,86.6),(0,1911,84.6),(0,1912,78.8),(0,1915,76.2),(0,1920,73.6),(0,1923,72.8),(0,1924,72.2),(0,1926,70),(0,1929,69.4),(0,1930,68),(0,1931,66.6),(0,1933,66),(0,1934,65.4),(0,1936,64.6),(0,1956,62),(0,1958,61.2),(0,1960,60.2),(0,1962,59.5),(0,1964,58.9),(0,1972,58.5),(0,1973,57.54),(0,1974,56.96),(0,1976,55.65),(0,1978,55.41),(0,1980,54.79),(0,1986,54.73),(0,1992,54.48),(0,1994,54.01),(0,2000,53.77),(0,2004,53.52)] 
       
M = find_fit(Data,model); M 
       
[a == 555.7167834308991, b == -0.2514636815413249, c ==
-9.797961450991004]
[a == 555.7167834308991, b == -0.2514636815413249, c == -9.797961450991004]
model = model.subs(a=M[0].rhs(),b=M[1].rhs(),c=M[2].rhs()); model 
       
(gender, year) |--> -9.797961450991004*gender -
0.2514636815413249*year + 555.7167834308991
(gender, year) |--> -9.797961450991004*gender - 0.2514636815413249*year + 555.7167834308991
points([d[1:] for d in Data if d[0]==1],figsize=5)+points([d[1:] for d in Data if d[0]==0],color='red')+plot(model(gender=1),(year,1900,2010),color='blue',legend_label="modeling men's times" )+plot(model(gender=0),(year,1900,2010),color='red',legend_label="modeling women's times") 
       

And this is exactly the same thing as solving a least squares problem!

A = matrix(RDF,[d[:-1] for d in Data if d[0]==1]); A 
       
31 x 2 dense matrix over Real Double Field (use the '.str()' method to
see the entries)
31 x 2 dense matrix over Real Double Field (use the '.str()' method to see the entries)
y = vector(matrix(RDF,[d[2:] for d in Data if d[0]==1]).transpose().list()); y 
       
(65.8, 65.6, 62.8, 61.6, 61.4, 60.4, 58.6, 57.4, 56.8, 56.6, 56.4, 55.9,
55.8, 55.4, 54.8, 54.6, 53.6, 52.9, 52.6, 52.2, 51.9, 51.22, 50.59,
49.44, 49.36, 49.24, 48.74, 48.42, 48.21, 48.18, 47.84)
(65.8, 65.6, 62.8, 61.6, 61.4, 60.4, 58.6, 57.4, 56.8, 56.6, 56.4, 55.9, 55.8, 55.4, 54.8, 54.6, 53.6, 52.9, 52.6, 52.2, 51.9, 51.22, 50.59, 49.44, 49.36, 49.24, 48.74, 48.42, 48.21, 48.18, 47.84)

Note that with "random" data, the inverse of $A^TA$ will exist.

(A.transpose()*A)^-1*A.transpose()*y 
       
(394.8373768074424, -0.17412930023829695)
(394.8373768074424, -0.17412930023829695)
plot(394-.174*x,(x,1900,2000),figsize=4)+points([d[1:] for d in Data if d[0]==1])