<jupyter_start><jupyter_text>You can read an overview of this Numerical Linear Algebra course in [this blog post](http://www.fast.ai/2017/07/17/num-lin-alg/).  The course was originally taught in the [University of San Francisco MS in Analytics](https://www.usfca.edu/arts-sciences/graduate-programs/analytics) graduate program.  Course lecture videos are [available on YouTube](https://www.youtube.com/playlist?list=PLtmWHNX-gukIc92m1K0P6bIOnZb-mg0hY) (note that the notebook numbers and video numbers do not line up, since some notebooks took longer than 1 video to cover).You can ask questions about the course on [our fast.ai forums](http://forums.fast.ai/c/lin-alg).  7. PageRank with Eigen Decompositions  Two Handy Tricks Here are two tools that we'll be using today, which are useful in general. 1\. [Psutil](https://github.com/giampaolo/psutil) is a great way to check on your memory usage.  This will be useful here since we are using a larger data set.<jupyter_code>import psutil
process = psutil.Process(os.getpid())
t = process.memory_info()
t.vms, t.rss
def mem_usage():
    process = psutil.Process(os.getpid())
    return process.memory_info().rss / psutil.virtual_memory().total
mem_usage()<jupyter_output><empty_output><jupyter_text>2\. [TQDM](https://github.com/tqdm/tqdm) gives you progress bars.<jupyter_code>from time import sleep
# Without TQDM
s = 0
for i in range(10):
    s += i
    sleep(0.2)
print(s)
# With TQDM
from tqdm import tqdm

s = 0
for i in tqdm(range(10)):
    s += i
    sleep(0.2)
print(s)<jupyter_output>100%|██████████| 10/10 [00:02<00:00,  4.96it/s]<jupyter_text>Motivation **Review**- What is SVD?- What are some applications of SVD?  Additional SVD Application An interesting use of SVD that I recently came across was as a step in the de-biasing of Word2Vec word embeddings, from [Quantifying and Reducing Stereotypes in Word Embeddings](https://arxiv.org/pdf/1606.06121.pdf)(Bolukbasi, et al).Word2Vec is a useful library released by Google that represents words as vectors.  The similarity of vectors captures semantic meaning, and analogies can be found, such as *Paris:France :: Tokyo: Japan*.  (source: [Vector Representations of Words](https://www.tensorflow.org/versions/r0.10/tutorials/word2vec/))However, these embeddings can implicitly encode bias, such as *father:doctor :: mother:nurse* and *man:computer programmer :: woman:homemaker*.One approach for de-biasing the space involves using SVD to reduce the dimensionality ([Bolukbasi paper](https://arxiv.org/pdf/1606.06121.pdf)).You can read more about bias in word embeddings:- [How Vector Space Mathematics Reveals the Hidden Sexism in Language](https://www.technologyreview.com/s/602025/how-vector-space-mathematics-reveals-the-hidden-sexism-in-language/)(MIT Tech Review)- [ConceptNet: better, less-stereotyped word vectors](https://blog.conceptnet.io/2017/04/24/conceptnet-numberbatch-17-04-better-less-stereotyped-word-vectors/)- [Semantics derived automatically from language corpora necessarily contain human biases](https://www.princeton.edu/~aylinc/papers/caliskan-islam_semantics.pdf) (excellent and very interesting paper!)  Ways to think about SVD - Data compression- SVD trades a large number of features for a smaller set of better features- All matrices are diagonal (if you use change of bases on the domain and range)  Perspectives on SVD We usually talk about SVD in terms of matrices, $$A = U \Sigma V^T$$ but we can also think about it in terms of vectors.  SVD gives us sets of orthonormal vectors ${v_j}$ and ${u_j}$ such that $$ A v_j = \sigma_j u_j $$$\sigma_j$ are scalars, called singular values**Q**: Does this remind you of anything?  Answer **Relationship between SVD and Eigen Decomposition**: the left-singular vectors of A are the eigenvectors of $AA^T$. The right-singular vectors of A are the eigenvectors of $A^T A$. The non-zero singular values of A are the square roots of the eigenvalues of $A^T A$ (and $A A^T$).SVD is a generalization of eigen decomposition. Not all matrices have eigen values, but ALL matrices have singular values.Let's forget SVD for a bit and talk about how to find the eigenvalues of a symmetric positive definite matrix...  Further resources on SVD - [SVD: image compression and least squares](http://andrew.gibiansky.com/blog/mathematics/cool-linear-algebra-singular-value-decomposition/)- [Image Compression with SVD](http://nbviewer.jupyter.org/gist/frankcleary/4d2bd178708503b556b0)- [The Extraordinary SVD](https://sites.math.washington.edu/~morrow/498_13/svd_applied.pdf)- [A Singularly Valuable Decomposition: The SVD of a Matrix](https://sites.math.washington.edu/~morrow/498_13/svd.pdf)  Today: Eigen Decomposition The best classical methods for computing the SVD are variants on methods for computing eigenvalues.  In addition to their links to SVD, Eigen decompositions are useful on their own as well.  Here are a few practical applications of eigen decomposition:- [rapid matrix powers](http://www.onmyphd.com/?p=eigen.decompositionh2_why)- [nth Fibonacci number](http://mathproofs.blogspot.com/2005/04/nth-term-of-fibonacci-sequence.html)- Behavior of ODEs- Markov Chains (health care economics, Page Rank)- [Linear Discriminant Analysis on Iris dataset](http://sebastianraschka.com/Articles/2014_python_lda.html) Check out the 3 Blue 1 Brown videos on [Change of basis](https://www.youtube.com/watch?v=P2LTAUO1TdA) and [Eigenvalues and eigenvectors](https://www.youtube.com/watch?v=PFDu9oVAE-g) "Eigenvalues are a way to see into the heart of a matrix... All the difficulties of matrices are swept away" -Strang  **Vocab**: A **Hermitian** matrix is one that is equal to it's own conjugate transpose.  In the case of real-valued matrices (which is all we are considering in this course), **Hermitian** means the same as **Symmetric**.**Relevant Theorems:**- If A is symmetric, then eigenvalues of A are real and $A = Q \Lambda Q^T$- If A is triangular, then its eigenvalues are equal to its diagonal entries  DBpedia Dataset Let's start with the **Power Method**, which finds one eigenvector.  *What good is just one eigenvector?* you may be wondering.  This is actually the basis for PageRank (read [The $25,000,000,000 Eigenvector: the Linear Algebra Behind Google](http://www.rose-hulman.edu/~bryan/googleFinalVersionFixed.pdf) for more info)Instead of trying to rank the importance of all websites on the internet, we are going to use a dataset of Wikipedia links from [DBpedia](http://wiki.dbpedia.org/).  DBpedia provides structured Wikipedia data available in 125 languages.  "*The full DBpedia data set features 38 million labels and abstracts in 125 different languages, 25.2 million links to images and 29.8 million links to external web pages; 80.9 million links to Wikipedia categories, and 41.2 million links to [YAGO](http://www.mpi-inf.mpg.de/departments/databases-and-information-systems/research/yago-naga/yago/) categories*" --[about DBpedia](http://wiki.dbpedia.org/about)Today's lesson is inspired by this [SciKit Learn Example](http://scikit-learn.org/stable/auto_examples/applications/wikipedia_principal_eigenvector.htmlsphx-glr-auto-examples-applications-wikipedia-principal-eigenvector-py)  Imports<jupyter_code>import os, numpy as np, pickle
from bz2 import BZ2File
from datetime import datetime
from pprint import pprint
from time import time
from tqdm import tqdm_notebook
from scipy import sparse

from sklearn.decomposition import randomized_svd
from sklearn.externals.joblib import Memory
from urllib.request import urlopen<jupyter_output><empty_output><jupyter_text>Downloading the data The data we have are:- **redirects**: URLs that redirect to other URLs- **links**: which pages link to which other pages Note: this takes a while<jupyter_code>PATH = 'data/dbpedia/'
URL_BASE = 'http://downloads.dbpedia.org/3.5.1/en/'
filenames = ["redirects_en.nt.bz2", "page_links_en.nt.bz2"]

for filename in filenames:
    if not os.path.exists(PATH+filename):
        print("Downloading '%s', please wait..." % filename)
        open(PATH+filename, 'wb').write(urlopen(URL_BASE+filename).read())
redirects_filename = PATH+filenames[0]
page_links_filename = PATH+filenames[1]<jupyter_output><empty_output><jupyter_text>Graph Adjacency Matrix We will construct a graph **adjacency matrix**, of which pages point to which.(source: [PageRank and HyperLink Induced Topic Search](https://www.slideshare.net/priyabrata232/page-rank-and-hyperlink))(source: [PageRank and HyperLink Induced Topic Search](https://www.slideshare.net/priyabrata232/page-rank-and-hyperlink))The power $A^2$ will give you how many ways there are to get from one page to another in 2 steps.  You can see a more detailed example, as applied to airline travel, [worked out in these notes](http://www.utdallas.edu/~jwz120030/Teaching/PastCoursesUMBC/M221HS06/ProjectFiles/Adjacency.pdf). We want to keep track of which pages point to which pages.  We will store this in a  square matrix, with a $1$ in position $(r, c)$ indicating that the topic in row $r$ points to the topic in column $c$You can read [more about graphs here](http://www.geeksforgeeks.org/graph-and-its-representations/).  Data Format One line of the file looks like:- `   .`In the below slice, the plus 1, -1 are to remove the<jupyter_code>DBPEDIA_RESOURCE_PREFIX_LEN = len("http://dbpedia.org/resource/")
SLICE = slice(DBPEDIA_RESOURCE_PREFIX_LEN + 1, -1)
def get_lines(filename): return (line.split() for line in BZ2File(filename))<jupyter_output><empty_output><jupyter_text>Loop through redirections and create dictionary of source to final destination<jupyter_code>def get_redirect(targ, redirects):
    seen = set()
    while True:
        transitive_targ = targ
        targ = redirects.get(targ)
        if targ is None or targ in seen: break
        seen.add(targ)
    return transitive_targ
def get_redirects(redirects_filename):
    redirects={}
    lines = get_lines(redirects_filename)
    return {src[SLICE]:get_redirect(targ[SLICE], redirects) 
                for src,_,targ,_ in tqdm_notebook(lines, leave=False)}
redirects = get_redirects(redirects_filename)
mem_usage()
def add_item(lst, redirects, index_map, item):
    k = item[SLICE]
    lst.append(index_map.setdefault(redirects.get(k, k), len(index_map)))
limit=119077682 #5000000
# Computing the integer index map
index_map = dict() # links->IDs
lines = get_lines(page_links_filename)
source, destination, data = [],[],[]
for l, split in tqdm_notebook(enumerate(lines), total=limit):
    if l >= limit: break
    add_item(source, redirects, index_map, split[0])
    add_item(destination, redirects, index_map, split[2])
    data.append(1)
n=len(data); n<jupyter_output><empty_output><jupyter_text>Looking at our data The below steps are just to illustrate what info is in our data and how it is structured.  They are not efficient. Let's see what type of items are in index_map:<jupyter_code>index_map.popitem()<jupyter_output><empty_output><jupyter_text>Let's look at one item in our index map: 1940_Cincinnati_Reds_Team_Issue has index $9991173$.  This only shows up once in the destination list:<jupyter_code>[i for i,x in enumerate(source) if x == 9991173]
source[119077649], destination[119077649]<jupyter_output><empty_output><jupyter_text>Now, we want to check which page is the source (has index $9991050$).  Note: usually you should not access a dictionary by searching for its values.  This is inefficient and not how dictionaries are intended to be used.<jupyter_code>for page_name, index in index_map.items():
    if index == 9991050:
        print(page_name)<jupyter_output>b'W711-2'<jupyter_text>We can see on Wikipedia that the Cincinati Red Teams Issue has [redirected to W711-2](https://en.wikipedia.org/wiki/W711-2):<jupyter_code>test_inds = [i for i,x in enumerate(source) if x == 9991050]
len(test_inds)
test_inds[:5]
test_dests = [destination[i] for i in test_inds]<jupyter_output><empty_output><jupyter_text>Now, we want to check which page is the source (has index 9991174):<jupyter_code>for page_name, index in index_map.items():
    if index in test_dests:
        print(page_name)<jupyter_output>b'Baseball'
b'Ohio'
b'Cincinnati'
b'Flash_Thompson'
b'1940'
b'1938'
b'Lonny_Frey'
b'Cincinnati_Reds'
b'Ernie_Lombardi'
b'Baseball_card'
b'James_Wilson'
b'Trading_card'
b'Detroit_Tigers'
b'Baseball_culture'
b'Frank_McCormick'
b'Bucky_Walters'
b'1940_World_Series'
b'Billy_Werber'
b'Ival_Goodman'
b'Harry_Craft'
b'Paul_Derringer'
b'Johnny_Vander_Meer'
b'Cigarette_card'
b'Eddie_Joost'
b'Myron_McCormick'
b'Beckett_Media'
b'Icarus_affair'
b'Ephemera'
b'Sports_card'
b'James_Turner'
b'Jimmy_Ripple'
b'Lewis_Riggs'
b'The_American_Card_Catalog'
b'Rookie_card'
b'Willard_Hershberger'
b'Elmer_Riddle'
b'Joseph_Beggs'
b'Witt_Guise'
b'Milburn_Shoffner'<jupyter_text>We can see that the items in the list appear in the wikipedia page:(Source: [Wikipedia](https://en.wikipedia.org/wiki/W711-2))  Create Matrix Below we create a sparse matrix using Scipy's COO format, and that convert it to CSR.**Questions**: What are COO and CSR?  Why would we create it with COO and then convert it right away?<jupyter_code>X = sparse.coo_matrix((data, (destination,source)), shape=(n,n), dtype=np.float32)
X = X.tocsr()
del(data,destination, source)
X
names = {i: name for name, i in index_map.items()}
mem_usage()<jupyter_output><empty_output><jupyter_text>Save matrix so we don't have to recompute<jupyter_code>pickle.dump(X, open(PATH+'X.pkl', 'wb'))
pickle.dump(index_map, open(PATH+'index_map.pkl', 'wb'))
X = pickle.load(open(PATH+'X.pkl', 'rb'))
index_map = pickle.load(open(PATH+'index_map.pkl', 'rb'))
names = {i: name for name, i in index_map.items()}
X<jupyter_output><empty_output><jupyter_text>Power method  Motivation An $n \times n$ matrix $A$ is **diagonalizable** if it has $n$ linearly independent eigenvectors $v_1,\, \ldots v_n$.Then any $w$ can be expressed $w = \sum_{j=1}^n c_j v_j $, for some scalars $c_j$.**Exercise:** Show that $$ A^k w = \sum_{j=1}^n c_j \lambda_j^k v_j$$**Question**: How will this behave for large $k$? This is inspiration for the **power method**.  Code<jupyter_code>def show_ex(v):
    print(', '.join(names[i].decode() for i in np.abs(v.squeeze()).argsort()[-1:-10:-1]))
?np.squeeze<jupyter_output><empty_output><jupyter_text>How to normalize a sparse matrix:<jupyter_code>S = sparse.csr_matrix(np.array([[1,2],[3,4]]))
Sr = S.sum(axis=0).A1; Sr<jupyter_output><empty_output><jupyter_text>[numpy.matrix.A1](https://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.A1.htmlnumpy.matrix.A1)<jupyter_code>S.indices
S.data
S.data / np.take(Sr, S.indices)
def power_method(A, max_iter=100):
    n = A.shape[1]
    A = np.copy(A)
    A.data /= np.take(A.sum(axis=0).A1, A.indices)

    scores = np.ones(n, dtype=np.float32) * np.sqrt(A.sum()/(n*n)) # initial guess
    for i in range(max_iter):
        scores = A @ scores
        nrm = np.linalg.norm(scores)
        scores /= nrm
        print(nrm)

    return scores<jupyter_output><empty_output><jupyter_text>**Question**: Why normalize the scores on each iteration?<jupyter_code>scores = power_method(X, max_iter=10)
show_ex(scores)
mem_usage()<jupyter_output><empty_output><jupyter_text>Comments Many advanced eigenvalue algorithms that are used in practice are variations on the power method.In Lesson 3: Background Removal, we used Facebook's [fast randomized pca/svd library, fbpca](https://github.com/facebook/fbpca).  Check out [the source code](https://github.com/facebook/fbpca/blob/master/fbpca.pyL1549) for the pca method we used. It uses the power method! **Further Study**- Check out [Google Page Rank, Power Iteration and the Second EigenValue of the Google Matrix](http://rstudio-pubs-static.s3.amazonaws.com/239261_8a607707294341c4b7e26acf728c28bd.html) for animations of the distribution as it converges.- The convergence rate of the power method is the ratio of the largest eigenvalue to the 2nd largest eigenvalue.  It can be speeded up by adding *shifts*.  To find eigenvalues other than the largest, a method called *deflation* can be used.  See Chapter 12.1 of Greenbaum & Chartier for more details.- *Krylov Subspaces*: In the Power Iteration, notice that we multiply by our matrix A each time, effectively calculating $$Ab,\,A^2b,\,A^3b,\,A^4b\, \ldots$$    The matrix with those vectors as columns is called the *Krylov matrix*, and the space spanned by those vectors is the *Krylov subspace*.  Keep this in mind for later.  Compare to SVD<jupyter_code>%time U, s, V = randomized_svd(X, 3, n_iter=3)
mem_usage()
# Top wikipedia pages according to principal singular vectors
show_ex(U.T[0])
show_ex(U.T[1])
show_ex(V[0])
show_ex(V[1])<jupyter_output>Democratic_Party_%28United_States%29, Republican_Party_%28United_States%29, Democratic-Republican_Party, United_States, Whig_Party_%28United_States%29, Federalist_Party, National_Republican_Party, Catholic_Church, Jacksonian_democracy<jupyter_text>**Exercise:** Normalize the data in various ways.  Don't overwrite the adjacency matrix, but instead create a new one.  See how your results differ. **Eigen Decomposition vs SVD:**- SVD involves 2 bases, eigen decomposition involves 1 basis- SVD bases are orthonormal, eigen basis generally not orthogonal- All matrices have an SVD, not all matrices (not even all square) have an eigen decomposition  QR Algorithm We used the power method to find the eigenvector corresponding to the largest eigenvalue of our matrix of Wikipedia links.  This eigenvector gave us the relative importance of each Wikipedia page (like a simplified PageRank).Next, let's look at a method for finding all eigenvalues of a symmetric, positive definite matrix.  This method includes 2 fundamental algorithms in numerical linear algebra, and is a basis for many more complex methods.[The Second Eigenvalue of the Google Matrix](https://nlp.stanford.edu/pubs/secondeigenvalue.pdf): has "implications for the convergence rate of the standard PageRank algorithm as the web scales, for the stability of PageRank to perturbations to the link structure of the web, for the detection of Google spammers, and for the design of algorithms to speed up PageRank".  Avoiding Confusion: QR Algorithm vs QR Decomposition The **QR algorithm** uses something called the **QR decomposition**.  Both are important, so don't get them confused.  The **QR decomposition** decomposes a matrix $A = QR$ into a set of orthonormal columns $Q$ and a triangular matrix $R$.  We will look at several ways to calculate the QR decomposition in a future lesson.  For now, just know that it is giving us an orthogonal matrix and a triangular matrix.  Linear Algebra Two matrices $A$ and $B$ are **similar** if there exists a non-singular matrix $X$ such that $$B = X^{-1}AX$$Watch this: [Change of Basis](https://www.youtube.com/watch?v=P2LTAUO1TdA&index=13&list=PLZHQObOWTQDPD3MizzM2xVFitgF8hE_ab)**Theorem**: If $X$ is non-singular, then $A$ and $X^{-1}AX$ have the same eigenvalues.  More Linear Algebra A **Schur factorization** of a matrix $A$ is a factorization:$$ A = Q T Q^*$$where $Q$ is unitary and $T$ is upper-triangular.**Question**: What can you say about the eigenvalues of $A$?**Theorem:** Every square matrix has a Schur factorization.  Other resources Review: [Linear combinations, span, and basis vectors](https://www.youtube.com/watch?v=k7RM-ot2NWY&index=3&list=PLZHQObOWTQDPD3MizzM2xVFitgF8hE_ab)See Lecture 24 for proofs of above theorems (and more!)  Algorithm The most basic version of the QR algorithm:    for k=1,2,...        Q, R = A        A = R @ Q        Under suitable assumptions, this algorithm converges to the Schur form of A!  Why it works Written again, only with subscripts:$A_0 = A$for $k=1,2,\ldots$   $\quad Q_k$, $R_k$ = $A_{k-1}$       $\quad A_k$ = $R_k Q_k$        We can think of this as constructing sequences of $A_k$, $Q_k$, and $R_k$.$$ A_k = Q_k \, R_k $$$$ Q_k^{-1} \, A_k = R_k$$Thus, $$ R_k Q_k = Q_k^{-1} \, A_k \, Q_k $$$$A_k = Q_k^{-1} Q_2^{-1} Q_1^{-1} A Q_1 Q_2 \dots Q_k$$Trefethen proves the following on page 216-217:$$A^k = Q_1 Q_2 \dots Q_k R_k R_{k-1}\dots R_1$$**Key**: The QR algorithm constructs orthonormal bases for successive powers $A^k$.  And remember the close relationship between powers of A and the eigen decomposition.To learn more, read up on *Rayleigh quotients*.  Pure QR<jupyter_code>n = 6
A = np.random.rand(n,n)
AT = A @ A.T
def pure_qr(A, max_iter=50000):
    Ak = np.copy(A)
    n = A.shape[0]
    QQ = np.eye(n)
    for k in range(max_iter):
        Q, R = np.linalg.qr(Ak)
        Ak = R @ Q
        QQ = QQ @ Q
        if k % 100 == 0:
            print(Ak)
            print("\n")
    return Ak, QQ<jupyter_output><empty_output><jupyter_text>Pure QR<jupyter_code>Ak, Q = pure_qr(A)<jupyter_output>[[ 2.65646  0.21386  0.16765  0.75256  0.61251  0.93518]
 [ 0.52744  0.47579  0.17052 -0.41086 -0.21182 -0.01876]
 [ 0.29923  0.06964  0.11173  0.1879  -0.29101  0.60032]
 [ 0.2274   0.46162 -0.26654  0.08899  0.24631  0.26447]
 [-0.06093  0.02892  0.34162  0.07533  0.02393 -0.05456]
 [-0.06025  0.02694 -0.11675 -0.00927 -0.11939 -0.00767]]


[[ 2.78023  0.52642  0.0395  -0.11135  0.1569   1.15184]
 [ 0.       0.18624 -0.297   -0.07256 -0.04537  0.27907]
 [ 0.       0.69328  0.34105 -0.12198  0.11029  0.0992 ]
 [-0.      -0.0494  -0.02057  0.09461  0.59466  0.09115]
 [-0.       0.00008 -0.02659 -0.40372  0.06542  0.38612]
 [-0.       0.       0.       0.      -0.      -0.11832]]


[[ 2.78023 -0.12185 -0.51401  0.17625 -0.07467  1.15184]
 [ 0.       0.2117  -0.70351  0.09974 -0.02986  0.00172]
 [ 0.       0.28284  0.32635 -0.0847  -0.08488 -0.29104]
 [-0.      -0.00068 -0.00088 -0.01282  0.54836  0.13447]
 [-0.       0.      -0.00102 -0.45718  0.16208 -0.37726]
 [-0.       0.       0.  [...]<jupyter_text>Let's compare to the eigenvalues:<jupyter_code>np.linalg.eigvals(A)<jupyter_output><empty_output><jupyter_text>Check that Q is orthogonal:<jupyter_code>np.allclose(np.eye(n), Q @ Q.T), np.allclose(np.eye(n), Q.T @ Q)<jupyter_output><empty_output><jupyter_text>This is really really slow.  Practical QR (QR with shifts) **Idea**: Instead of factoring $A_k$ as $Q_k R_k$, 1. Get the QR factorization $$A_k - s_k I = Q_k R_k$$2. Set $$A_{k+1} = R_k Q_k + s_k I$$Choose $s_k$ to approximate an eigenvalue of $A$.  We'll use $s_k = A_k(m,m)$. The idea of adding shifts to speed up convergence shows up in many algorithms in numerical linear algebra (including the power method, inverse iteration, and Rayleigh quotient iteration).     Homework: Add shifts to the QR algorithm<jupyter_code>#Exercise: Add shifts to the QR algorithm
#Exercise: def practical_qr(A, iters=10):
#Exercise:     return Ak, Q<jupyter_output><empty_output><jupyter_text>Practical QR<jupyter_code>Ak, Q = practical_qr(A, 10)<jupyter_output>[ 5.16264  2.53603  0.31923  0.35315  0.97569  0.43615]
[ 7.99381  0.05922  0.34478  0.29482  0.79026  0.29999]
[ 8.00491  0.04358  0.89735  0.26386  0.26182  0.31135]
[ 8.00493  0.13648  0.91881  0.14839  0.24313  0.33115]
[ 8.00493  0.43377  0.62809  0.13429  0.24592  0.33589]
[ 8.00493  0.81058  0.25128  0.13297  0.24722  0.3359 ]
[ 8.00493  0.98945  0.07221  0.13292  0.24747  0.3359 ]
[ 8.00493  1.0366   0.02497  0.13296  0.24751  0.3359 ]
[ 8.00493  1.04688  0.01465  0.13299  0.24752  0.3359 ]
[ 8.00493  1.04902  0.0125   0.13301  0.24753  0.3359 ]<jupyter_text>Check that Q is orthogonal:<jupyter_code>np.allclose(np.eye(n), Q @ Q.T), np.allclose(np.eye(n), Q.T @ Q)<jupyter_output><empty_output><jupyter_text>Let's compare to the eigenvalues:<jupyter_code>np.linalg.eigvals(A)<jupyter_output><empty_output><jupyter_text>**Problem**: This is better than the unshifted version (which wasn't even guaranteed to converge), but is still really slow!  In fact, it is $\mathcal{O}(n^4)$, which is awful. In the case of symmetric matrices, it's $\mathcal{O}(n^3)$ However, if you start with a **Hessenberg matrix** (zeros below the first subdiagonal), it's faster: $\mathcal{O}(n^3)$, and $\mathcal{O}(n^2)$ if symmetric.  A Two-Phase Approach In practice, a two phase approach is used to find eigenvalues:1. Reduce the matrix to *Hessenberg* form (zeros below the first subdiagonal)2. Iterative process that causes Hessenberg to converge to a *triangular* matrix.  The eigenvalues of a triangular matrix are the values on the diagonal, so we are finished!(source: Trefethen, Lecture 25)In the case of a Hermitian matrix, this approach is even faster, since the intermediate step is also Hermitian (and a Hermitian Hessenberg is *tridiagonal*).(source: Trefethen, Lecture 25)Phase 1 reaches an exact solution in a finite number of steps, whereas Phase 2 theoretically never reaches the exact solution.We've already done step 2: the QR algorithm.  Remember that it would be possible to just use the QR algorithm, but ridiculously slow.  Arnoldi Iteration We can use the Arnoldi iteration for phase 1 (and the QR algorithm for phase 2).  Initializations<jupyter_code>import numpy as np
n = 5
A0 = np.random.rand(n,n)  #.astype(np.float64)
A = A0 @ A0.T

np.set_printoptions(precision=5, suppress=True)<jupyter_output><empty_output><jupyter_text>Linear Algebra Review: Projections When vector $\mathbf{b}$ is projected onto a line $\mathbf{a}$, its projection $\mathbf{p}$ is the part of $\mathbf{b}$ along that line $\mathbf{a}$.Let's look at interactive graphic (3.4) for [section 3.2.2: Projections](http://immersivemath.com/ila/ch03_dotproduct/ch03.html) of the [Immersive Linear Algebra online book](http://immersivemath.com/ila/index.html).(source: [Immersive Math](http://immersivemath.com/ila/ch03_dotproduct/ch03.html))And here is what it looks like to project a vector onto a plane:(source: [The Linear Algebra View of Least-Squares Regression](https://medium.com/@andrew.chamberlain/the-linear-algebra-view-of-least-squares-regression-f67044b7f39b))When vector $\mathbf{b}$ is projected onto a line $\mathbf{a}$, its projection $\mathbf{p}$ is the part of $\mathbf{b}$ along that line $\mathbf{a}$.  So $\mathbf{p}$ is some multiple of $\mathbf{a}$. Let $\mathbf{p} = \hat{x}\mathbf{a}$ where $\hat{x}$ is a scalar.  Orthogonality **The key to projection is orthogonality:** The line *from* $\mathbf{b}$ to $\mathbf{p}$ (which can be written $\mathbf{b} - \hat{x}\mathbf{a}$) is perpendicular to $\mathbf{a}$.This means that $$ \mathbf{a} \cdot (\mathbf{b} -  \hat{x}\mathbf{a}) = 0 $$and so $$\hat{x} = \frac{\mathbf{a} \cdot \mathbf{b}}{\mathbf{a} \cdot \mathbf{a}} $$  The Algorithm **Motivation**:We want orthonormal columns in $Q$ and a Hessenberg $H$ such that $A Q = Q H$.Thinking about it iteratively, $$ A Q_n = Q_{n+1} H_n $$ where $Q_{n+1}$ is $n\times n+1$ and $H_n$ is $n+1 \times n$.  This creates a solvable recurrence relation.(source: Trefethen, Lecture 33) **Pseudo-code for Arnoldi Algorithm**    Start with an arbitrary vector (normalized to have norm 1) for first col of Q    for n=1,2,3...        v = A @ nth col of Q        for j=1,...n            project v onto q_j, and subtract the projection off of v            want to capture part of v that isn't already spanned by prev columns of Q            store coefficients in H        normalize v, and then make it the (n+1)th column of Q Notice that we are multiplying A by the previous vector in Q and removing the components that are not orthogonal to the existing columns of Q.**Question:** Repeated multiplications of A?  Does this remind you of anything?  Answer:<jupyter_code>#Exercise Answer
The *Power Method* involved iterative multiplications by A as well!<jupyter_output><empty_output><jupyter_text>About how the Arnoldi Iteration works - With the Arnoldi Iteration, we are finding an orthonormal basis for the *Krylov subspace*. The Krylov matrix $$ K = \left[b \; Ab \; A^2b \; \dots \; A^{n-1}b \right]$$has a QR factorization$$K = QR$$and that is the same $Q$ that is being found in the Arnoldi Iteration.  Note that the Arnoldi Iteration does not explicity calculate $K$ or $R$.- Intuition: K contains good information about the largest eigenvalues of A, and the QR factorization reveals this information by peeling off one approximate eigenvector at a time. The Arnoldi Iteration is two things:1. the basis of many of the iterative algorithms of numerical linear algebra2. a technique for finding eigenvalues of nonhermitian matrices(Trefethen, page 257)**How Arnoldi Locates Eigenvalues**1. Carry out Arnoldi iteration2. Periodically calculate the eigenvalues (called *Arnoldi estimates* or *Ritz values*) of the Hessenberg H, using the QR algorithm3. Check at whether these values are converging.  If they are, they're probably eigenvalues of A.  Implementation<jupyter_code># Decompose square matrix A @ Q ~= Q @ H
def arnoldi(A):
    m, n = A.shape
    assert(n <= m)
    
    # Hessenberg matrix
    H = np.zeros([n+1,n]) #, dtype=np.float64)
    # Orthonormal columns
    Q = np.zeros([m,n+1]) #, dtype=np.float64)
    # 1st col of Q is a random column with unit norm
    b = np.random.rand(m)
    Q[:,0] = b / np.linalg.norm(b)
    for j in range(n):
        v = A @ Q[:,j]
        for i in range(j+1):
            #This comes from the formula for projection of v onto q.
            #Since columns q are orthonormal, q dot q = 1
            H[i,j] = np.dot(Q[:,i], v)
            v = v - (H[i,j] * Q[:,i])
        H[j+1,j] = np.linalg.norm(v)
        Q[:,j+1] = v / H[j+1,j]
        
        # printing this to see convergence, would be slow to use in practice
        print(np.linalg.norm(A @ Q[:,:-1] - Q @ H))
    return Q[:,:-1], H[:-1,:]
Q, H = arnoldi(A)<jupyter_output>8.59112969133
4.45398729097
0.935693639985
3.36613943339
0.817740180293<jupyter_text>Check that H is tri-diagonal:<jupyter_code>H<jupyter_output><empty_output><jupyter_text>Exercise Write code to confirm that:1. AQ = QH2. Q is orthonormal  Answer:<jupyter_code>#Exercise:
np.allclose(A @ Q, Q @ H)
#Exercise:
np.allclose(np.eye(len(Q)), Q.T @ Q)<jupyter_output><empty_output><jupyter_text>General Case: **General Matrix**: Now we can do this on our general matrix A (not symmetric).  In this case, we are getting a Hessenberg instead of a Tri-diagonal<jupyter_code>Q0, H0 = arnoldi(A0)<jupyter_output>1.44287067354
1.06234006889
0.689291414367
0.918098818651
4.7124490411e-16<jupyter_text>Check that H is Hessenberg:<jupyter_code>H0
np.allclose(A0 @ Q0, Q0 @ H0)
np.allclose(np.eye(len(Q0)), Q0.T @ Q0), np.allclose(np.eye(len(Q0)), Q0 @ Q0.T)<jupyter_output><empty_output><jupyter_text>Putting it all together<jupyter_code>def eigen(A, max_iter=20):
    Q, H = arnoldi(A)
    Ak, QQ = practical_qr(H, max_iter)
    U = Q @ QQ
    D = np.diag(Ak)
    return U, D
n = 10
A0 = np.random.rand(n,n)
A = A0 @ A0.T
U, D = eigen(A, 40)
D
np.linalg.eigvals(A)
np.linalg.norm(U @ np.diag(D) @ U.T - A)
np.allclose(U @ np.diag(D) @ U.T, A, atol=1e-3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1149458811329315, 'fraction_numerical': 0.051137468403655455, 'mean_word_length': 4.493857931280043, 'pattern_counts': {'":': 0, '<': 131, '<?xml version=': 0, '>': 131, 'https://': 25, 'lorem ipsum': 0, 'www.': 17, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '282035', 'license': 'Apache-2.0', 'n_tokens_mistral': 11609, 'n_tokens_neox': 10373, 'n_words': 3438, 'path': 'numerical-linear-algebra/nbs/7. PageRank with Eigen Decompositions.ipynb', 'repo_name': 'nimacgit/courses', 'size': '30858'}