Be sure to make a copy!!!!
This notebook walks you through the operations required to compute a low rank approximation of a matrix $ B $. We will create a matrix $A$ whose column space will be used in the approximation of $B$.
We start by creating a random $ m \times n $ matrix $ B $. We then take $ k $ columns of $ B $ to be matrix $ A $, whose columns will be used in the approximation $ B \approx A V $. (In the text and the videos, we talk about computing $ W $ so that $ B \approx A W^T $. Here we find it more convenient to compute the transpose of that matrix instead. We call it $ V $ to distinguish it from $ W $. So, $ W = V^T $.)
$ V $ is computed as $ ( A^T A )^{-1} A^T B $.
import numpy as np
import laff
import flame
m = 8
n = 8
k = 3
# Random matrix of size mxn
B = np.matrix( np.random.rand( m, n ) )
# A is k columns of B taken at even intervals
if 2*k <= n: #k is less than half of n
interval = np.ceil( n/k )
A = B[ :, ::interval ] # These are slices in python.
# This says take all rows of B, and columns
# from 0 to the end at interval steps
else:
A = B[ :, :k] #If k is greater than half of n, then just take the first k columns
print( 'A = ' )
print( A )
print( 'B = ' )
print( B )
We start the process of computing $( A^T A )^{-1} A^T B$ by computing $ A^T A $ and storing the result in a matrix, $C$. In this implementation, we ignore symmetry.
C = np.transpose( A ) * A
print( 'C = ' )
print( C )
Now, form $ V = A^T B $, notice that we are not done forming $V$ after this step.
V = np.transpose( A ) * B
Instead of computing $ C^{-1} = ( A^T A )^{-1} $ explicitly, we notice that we can instead store the $ L $ and $ U $ factorization of $C$ in $ C $ and then just solve $ L ( U X ) = V $. First we will overwrite $ V $ with the result of solving $ L Z = V $, and then we will overwrite $ V $ with the result of solving $ U X = V $.
Copy your LU_unb_var5
routine from Notebook 6.3: Solving A x b via LU Factorization and Triangular Solves
import flame
import laff
def LU_unb_var5(A):
ATL, ATR, \
ABL, ABR = flame.part_2x2(A, \
0, 0, 'TL')
while ATL.shape[0] < A.shape[0]:
A00, a01, A02, \
a10t, alpha11, a12t, \
A20, a21, A22 = flame.repart_2x2_to_3x3(ATL, ATR, \
ABL, ABR, \
1, 1, 'BR')
#------------------------------------------------------------#
laff.invscal( alpha11, a21 ) # a21 := a21 / alpha11
laff.ger( -1.0, a21, a12t, A22 ) # A22 := A22 - a21 * a12t
#------------------------------------------------------------#
ATL, ATR, \
ABL, ABR = flame.cont_with_3x3_to_2x2(A00, a01, A02, \
a10t, alpha11, a12t, \
A20, a21, A22, \
'TL')
flame.merge_2x2(ATL, ATR, \
ABL, ABR, A)
Now run LU_unb_var5
on the matrix $C$ to store $L$ and $U$ in it.
LU_unb_var5( C )
Solve $L ( U X ) = V$, overwriting $V$ where $U$ and $L$ are stored in the upper and the strictly lower portions of $C$ respectively.
laff.trsm('Lower triangular', 'No transpose', 'Unit diagonal', C, V)
laff.trsm('Upper triangular', 'No transpose', 'Nonunit diagonal', C, V)
The $ j $th column of $ A V $ now equals the projection of the $ j $th column of $ B $ onto the column space of $ A $, $ {\cal C}( A ) $.
A couple of notes:
The matrix $ A^T A $ is symmetric positive definite. As a result, one does not need to pivot when performing LU factorization. (The reason for this is beyond the scope of this course.)
One could use what is called a "Symmetric rank-k update" operation to compute only the lower (or upper) triangular part of $ A^T A $. This would (approximately) halve the number of floating point operations that are required.
In one of the enrichments, 8.4.2, we discussed the Cholesky factorization of a symmetric positive definite matrix. One should ideally use that here, since it also takes advantage of symmetry.
This would then leave us with $ L $, a lower triangular matrix, such that $ C = A^T A = L L^T $. Computing $ V $ would then require the steps
The above computation should be implemented as the routine RankKApprox( B, k )
where $ B $ is the $ m \times n $ matrix to be approximated, and $k$ is the rank of the eventual approximation that will be returned by the method.
def RankKApprox( B, k ):
m,n = B.shape # How many rows and columns does B have?
# A is k columns of B taken at even intervals
if 2*k <= n: #k is less than half of n
interval = np.ceil( n/k )
A = B[ :, ::interval ] # These are slices in python.
# This says take all rows of B, and columns
# from 0 to the end at interval steps
else:
A = B[ :, :k] #If k is greater than half of n, then just take the first k columns
# C = A^T A
C = np.transpose( A ) * A
# W = A^T B
W = np.transpose( A ) * B
# Overwrite C with its LU factorization
LU_unb_var5( C )
# Solve L(UX) = W, overwriting W with X
laff.trsm('Lower triangular', 'No transpose', 'Unit diagonal', C, W)
laff.trsm('Upper triangular', 'No transpose', 'Nonunit diagonal', C, W)
return A * W
Now that we have implemented routines to create low rank approximations to matrices we will explore what a rank k approximation to an image looks like. Each pixel in an image can be thought of as a value in a matrix. For a grayscale image, this value corresponds to how black or white it is on a relative scale.
We will use two techniques for these approximations. First, the normal approximation developed above and second, the SVD which is a very useful matrix decomposition that guarantees the best approximation given $k$ columns. The SVD might take a while to compute, so don't panic if one of the code blocks takes a bit to complete.
Try experimenting with the number of columns below by changing the numCols
variable.
If you want to use your own images, make sure that they are in the png
format and just place them in your notebooks directoy. Then change the filename
variable to reflect the file name of your image.
%pylab inline
from im_approx import *
# Try varying the number of columns used for the approximation
numCols = 40
# Load an image
filename = 'building.png'
img = np.matrix(read_image( filename ))
# Make the approximations using normal equations and the SVD
# This step might take a while as it is a lot of computation.
normalApprox, SVDApprox = create_approximations( img, k=numCols, approximator=RankKApprox )
#Now plot the approximations that we have created
# Note that we're having some issues with our
# approximations being somewhat darker than the
# real image and are investigating.
plot_approximations( img, normalApprox, SVDApprox, k=numCols )
# C = A^T A: This comment corresponds to computing $C = A^T A$. Try using numpy's built in transpose method by calling np.transpse(A)
.
# W = A^T B: This comment corresponds to computing $W = A^T B$. See above for a hint.
# Overwrite C with its LU factorization: Use your implementation of LU_unb_var5 from an earlier notebook for this.
# Solve L(UX) = W, overwriting W with X: Use laff.trsm
to do this. Recall that "trsm" means triangular solve with multiple right hand sides.