Is Q mxn or nxn? It is both! This is useful!¶

In [13]:

using LinearAlgebra
A = rand(4,2)
Q,R = qr(A)
Q

Out[13]:

4×4 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.36166    0.916348   -0.135025   -0.106193 
 -0.875369  -0.396235   -0.176529   -0.21346  
 -0.216      0.0556349   0.974526   -0.0234073
 -0.237226   0.0144527  -0.0300856   0.970881

It looks 4x4 !!¶

In [14]:

size(Q)

Out[14]:

(4, 4)

but wait! I can multiply by a vector of size 2¶

In [15]:

Q * [1,2]

Out[15]:

4-element Array{Float64,1}:
  1.4710364530489386 
 -1.6678387032016588 
 -0.10473038358204545
 -0.2083203318415276

what about size 3? !! (answer: no)¶

In [16]:

Q * [1,2,3]

DimensionMismatch("vector must have length either 4 or 2")

Stacktrace:
 [1] *(::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Int64,1}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/qr.jl:563
 [2] top-level scope at In[16]:1

what about size 4? (answer: yes)¶

In [20]:

Q * [1,2,3,4]

Out[20]:

4-element Array{Float64,1}:
  0.6411891565050274
 -3.051264992118778 
  2.7252182423170517
  3.5849468538868363

What's going on??¶

Q is not stored as elements, it is stored in a more compact form known as a WY representation which we do not cover in 18.06. This form not only saves memory, but allows us to complete the tall-skinny mxn Q into a full square orthogonal Q.

The "extra" vectors are an orthonormal set of vectors that are orthogonal to the column space of A. This is associated with the left nullspace of A.

In [1]:

A = rand(4,4)

Out[1]:

4×4 Array{Float64,2}:
 0.0486726  0.422254  0.0384858  0.324084
 0.700173   0.479445  0.570882   0.340941
 0.868803   0.791772  0.5067     0.294249
 0.550927   0.680756  0.133569   0.127481

In [2]:

A[:,2:end]

Out[2]:

4×3 Array{Float64,2}:
 0.422254  0.0384858  0.324084
 0.479445  0.570882   0.340941
 0.791772  0.5067     0.294249
 0.680756  0.133569   0.127481

In [3]:

A[:,2:4]

Out[3]:

4×3 Array{Float64,2}:
 0.422254  0.0384858  0.324084
 0.479445  0.570882   0.340941
 0.791772  0.5067     0.294249
 0.680756  0.133569   0.127481

In [4]:

A[:,[2,3,4]]

Out[4]:

4×3 Array{Float64,2}:
 0.422254  0.0384858  0.324084
 0.479445  0.570882   0.340941
 0.791772  0.5067     0.294249
 0.680756  0.133569   0.127481

Six Cases (kind of)¶

Square (rank=n rank<n)
Tall skinny (rank=n rank<n)
short wide (rank=m rank<m)

In [3]:

using LinearAlgebra
A = rand(3,3) # surely full rank
Q,R = qr(A)

Out[3]:

LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.688382   0.256761  -0.678383
 -0.536173  -0.810009   0.237495
 -0.488517   0.527218   0.695264
R factor:
3×3 Array{Float64,2}:
 -1.25521  -0.756215   -1.14882 
  0.0       0.0596418  -0.293094
  0.0       0.0        -0.393284

Q is square orthogonal , R is square, non-singular¶

In [6]:

v = A[:,1]

Out[6]:

3-element Array{Float64,1}:
 0.8640659111629196
 0.6730109100970267
 0.6131928969286227

In [7]:

norm(v)

Out[7]:

1.2552127759018195

In [8]:

v/norm(v)

Out[8]:

3-element Array{Float64,1}:
 0.6883820239497828
 0.5361727692848695
 0.4885170934370617

In [9]:

## what if it is square not full rank?

In [10]:

A = [ 1 2 3
      2 4 6
      1 1 1 ]

Out[10]:

3×3 Array{Int64,2}:
 1  2  3
 2  4  6
 1  1  1

In [11]:

rank(A)

Out[11]:

In [12]:

svdvals(A)

Out[12]:

3-element Array{Float64,1}:
 8.519782928662694    
 0.6428832308185806   
 7.528243337817958e-16

In [13]:

Q,R = qr(A)

Out[13]:

LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.408248  -0.182574  -0.894427   
 -0.816497  -0.365148   0.447214   
 -0.408248   0.912871  -8.88178e-16
R factor:
3×3 Array{Float64,2}:
 -2.44949  -4.49073   -6.53197    
  0.0      -0.912871  -1.82574    
  0.0       0.0        1.55431e-15

In [14]:

Q'Q

Out[14]:

3×3 Array{Float64,2}:
  1.0          -5.55112e-17  1.96064e-16
 -5.55112e-17   1.0          4.10453e-16
  1.96064e-16   4.10453e-16  1.0

if A is square and A = Q*R if A is singular, Q is still orthogonal, but R may not be invertible¶

(AB)⁻¹ = B⁻¹ * A⁻¹¶

Hence if the product of two matrices is not invertible then at least one of the two is not invertible¶

An orthogonal matrix has an inverse , i.e. Qᵀ¶

In [17]:

Q,R = qr(zeros(3,3))

Out[17]:

LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
R factor:
3×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

Tall Skinny¶

In [25]:

A = rand(6,3)

Out[25]:

6×3 Array{Float64,2}:
 0.253612  0.904386  0.54806 
 0.880715  0.814656  0.227191
 0.713426  0.799345  0.3588  
 0.502016  0.225784  0.623123
 0.64694   0.304476  0.987695
 0.481788  0.806257  0.262103

In [26]:

svdvals(A)

Out[26]:

3-element Array{Float64,1}:
 2.471516297197189 
 0.8480652908171992
 0.5015253422131328

In [27]:

rank(A)

Out[27]:

In [28]:

Q,R = qr(A);

In [29]:

Out[29]:

3×3 Array{Float64,2}:
 -1.50054  -1.47672   -1.11502  
  0.0       0.856476  -0.0306984
  0.0       0.0        0.817942

In [31]:

Q = Q[:,1:3]

Out[31]:

6×3 Array{Float64,2}:
 -0.169013   0.764528    0.468342
 -0.586931  -0.0608065  -0.52463 
 -0.475446   0.113538   -0.205207
 -0.334556  -0.313218    0.293994
 -0.431137  -0.387863    0.605252
 -0.321076   0.38777    -0.102697

In [32]:

Q*R ≈ A

Out[32]:

true

In [33]:

Q'Q

Out[33]:

3×3 Array{Float64,2}:
 1.0           3.46945e-16   9.71445e-17
 3.46945e-16   1.0          -2.35922e-16
 9.71445e-17  -2.35922e-16   1.0

In [34]:

Q*Q'

Out[34]:

6×6 Array{Float64,2}:
  0.832412   -0.192995   0.0710531  -0.0452299   0.0598008   0.30263  
 -0.192995    0.623423   0.379808    0.0611691  -0.040901    0.218749 
  0.0710531   0.379808   0.281049    0.0631715   0.0367432   0.217755 
 -0.0452299   0.0611691  0.0631715   0.296466    0.443666   -0.044231 
  0.0598008  -0.040901   0.0367432   0.443666    0.702647   -0.0741314
  0.30263     0.218749   0.217755   -0.044231   -0.0741314   0.264003

if A is tall-skinny, Q matches A (well julia has something here, but don't worry) and R is invertible square¶

What is this good for?¶

this is how least squares really happens

In [37]:

b = rand(6)

Out[37]:

6-element Array{Float64,1}:
 0.0478792078122805 
 0.27000278307530534
 0.31372385185711127
 0.8684428375752447 
 0.6468058271127577 
 0.4412365841937491

In [38]:

A\b

Out[38]:

3-element Array{Float64,1}:
  0.6267259740756708 
 -0.32726590196323513
  0.5108893852341747

In [41]:

R \( Q'b)

Out[41]:

3-element Array{Float64,1}:
  0.6267259740756709
 -0.3272659019632352
  0.5108893852341747

In [46]:

M = rand(4,4)

Out[46]:

4×4 Array{Float64,2}:
 0.0493362   0.574836  0.88775   0.470262
 0.00808782  0.711997  0.607727  0.476851
 0.917347    0.838914  0.104311  0.97497 
 0.723469    0.645267  0.240534  0.305889

In [53]:

A = [ 1 2 3;2 4 6; 3 6 9; 1 1 1; 1 1 1; 1 1 1]

Out[53]:

6×3 Array{Int64,2}:
 1  2  3
 2  4  6
 3  6  9
 1  1  1
 1  1  1
 1  1  1

In [54]:

rank(A)

Out[54]:

In [55]:

svdvals(A)

Out[55]:

3-element Array{Float64,1}:
 14.274567269631934  
  1.1120832993767427 
  3.2684080872602e-16

In [56]:

Q,R = qr(A)
Q = Q[:,1:3]

Out[56]:

6×3 Array{Float64,2}:
 -0.242536  -0.112272  -0.512343   
 -0.485071  -0.224544   0.783521   
 -0.727607  -0.336817  -0.351567   
 -0.242536   0.523937   1.11022e-16
 -0.242536   0.523937   1.11022e-16
 -0.242536   0.523937   1.11022e-16

In [57]:

Q*R ≈ A

Out[57]:

true

In [58]:

Out[58]:

6×3 Array{Float64,2}:
 -0.242536  -0.112272  -0.512343   
 -0.485071  -0.224544   0.783521   
 -0.727607  -0.336817  -0.351567   
 -0.242536   0.523937   1.11022e-16
 -0.242536   0.523937   1.11022e-16
 -0.242536   0.523937   1.11022e-16

In [59]:

Out[59]:

3×3 Array{Float64,2}:
 -4.12311  -7.5186   -10.9141     
  0.0      -1.57181   -3.14362    
  0.0       0.0       -4.00297e-16

In [60]:

Q'Q

Out[60]:

3×3 Array{Float64,2}:
  1.0          -5.55112e-17   8.57529e-17
 -5.55112e-17   1.0          -1.9783e-17 
  8.57529e-17  -1.9783e-17    1.0

summary of the four cases so far¶

for square and tall-skinny A = QR where the shape of Q is the shape of A and R is invertible or not, depending on whether A isfull rank or not

now for the short wide case¶

In [61]:

A = rand(3,6)

Out[61]:

3×6 Array{Float64,2}:
 0.864337  0.974244  0.508762  0.664655  0.556382  0.547616
 0.378073  0.491839  0.518706  0.902492  0.584839  0.389362
 0.555097  0.702286  0.426364  0.582923  0.553117  0.325505

In [62]:

rank(A)

Out[62]:

In [63]:

Q,R = qr(A);
R

Out[63]:

3×6 Array{Float64,2}:
 -1.0946  -1.29532    -0.797116  -1.13217   -0.921841  -0.731974 
  0.0     -0.0799796  -0.221097  -0.428477  -0.300855  -0.0683788
  0.0      0.0        -0.159382  -0.361468  -0.131273  -0.130301

In [64]:

Out[64]:

3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.789637   0.60756   -0.0856996
 -0.345398  -0.555592  -0.756318 
 -0.507123  -0.567616   0.648566

In [65]:

Q'Q

Out[65]:

3×3 Array{Float64,2}:
 1.0          2.22045e-16  1.66533e-16
 2.22045e-16  1.0          1.11022e-16
 1.66533e-16  1.11022e-16  1.0

A = Q * [R R₂] where R is non-singular upper triangular (full row rank)¶

In [66]:

A =[ 1 2 3;2 4 6; 3 6 9; 1 1 1; 1 1 1; 1 1 1]'

Out[66]:

3×6 Adjoint{Int64,Array{Int64,2}}:
 1  2  3  1  1  1
 2  4  6  1  1  1
 3  6  9  1  1  1

In [67]:

rank(A)

Out[67]:

In [68]:

Q,R = qr(A);

In [69]:

Out[69]:

3×6 Array{Float64,2}:
 -3.74166  -7.48331      -11.225        -1.60357   -1.60357   -1.60357 
  0.0       1.98603e-15    3.97205e-15   0.622821   0.622821   0.622821
  0.0       0.0           -3.9443e-31   -0.201656  -0.201656  -0.201656

In [70]:

Out[70]:

3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.267261   0.956183    0.119523
 -0.534522  -0.0439019  -0.844013
 -0.801784  -0.28946     0.522835

A = Q * [R R₂] where R is singular upper triangular (not full row rank)¶

For a matrix Q to satisfy Q'Q =I then Q must have no more columns than rows¶

example if such a Q existed for 3x6¶

then you would have 6 linearly independent vectors in R³

Qx=0 # ask is x=0 the only solution? Q'Qx=0=Ix so x=0 so the columns are linearly indpendent

In [ ]: