![\"N\times N\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-1f1ba6baee657f44c06c9a02c1d2fe7d_l3.png\" "\\"Rendered")
positive semidefinite (psd) matrix ![\"A\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-11f4f587954b361e7d78940f65b8d70d_l3.png\" "\\"Rendered")
. A special case was the column<\/em> Nystr\u00f6m approximation, defined to be1<\/a><\/sup>We use Matlab index notation<\/a> to indicate submatrices<\/a> of .<\/span>\n\n\n\n (Nys) <\/span> <\/span> where With the column Nystr\u00f6m approximation presented just as such, many questions remain:<\/p>\n\n\n\n In this post, we will answer all of these questions by drawing a connection between low-rank approximation by Nystr\u00f6m approximation and solving linear systems of equations by Gaussian elimination<\/a>. The connection between these two seemingly unrelated areas of matrix computations will pay dividends, leading to effective algorithms to compute Nystr\u00f6m approximations by the (partial) Cholesky factorization<\/a> of a positive (semi)definite matrix<\/a> and an elegant description of the residual of the Nystr\u00f6m approximation as the Schur complement<\/a>.<\/p>\n\n\n\n Suppose we want solve the system of linear equations For real symmetric positive definite matrix<\/a> The Cholesky factorization can be easily generalized to allow the matrix Here’s one way of computing the Cholesky factorization using recursion<\/em><\/a>. Write the matrix <\/span> <\/span> Our first step will be “block Cholesky factorize” the matrix The core idea of Gaussian elimination<\/a> is to combine rows of a matrix to introduce zero entries. For our case, observe that multiplying the first block row of <\/span> <\/span> Isolating <\/span> <\/span> yields the block triangular factorization<\/p>\n\n\n\n <\/span> <\/span> We’ve factored (1) <\/span> <\/span> This is a block version of the Cholesky decomposition of the matrix We’ve met the second of our main characters, the Schur complement<\/a><\/p>\n\n\n\n (Sch) <\/span> <\/span> This seemingly innocuous combination of matrices has a tendency to show up in surprising places when one works with matrices.4<\/a><\/sup>See my post on the Schur complement<\/a> for more examples.<\/span> It’s appearance in any one place is unremarkable, but the shear ubiquity of The Schur complement enjoys the following property:5<\/a><\/sup>This property is a consequence of equation (1) together with the conjugation rule for positive (semi)definiteness, which I discussed in this previous post<\/a>.<\/span>\n\n\n\n Positivity of the Schur complement:<\/strong> If As a consequence of this property, we conclude that both With the positive definiteness of the Schur complement in hand, we now return to our Cholesky factorization algorithm. Continue by recursively6<\/a><\/sup>As always with recursion, one needs to specify the base case. For us, the base case is just that Cholesky decomposition of a <\/span> <\/span> Inserting these into the block <\/span> <\/span> Voil\u00e0, we have obtained a Cholesky factorization of a positive definite matrix By unwinding the recursions (and always choosing the top left block <\/span> <\/span> <\/span> <\/span> This iterative algorithm is how Cholesky factorization is typically presented in textbooks<\/a>.<\/p>\n\n\n\n Our motivating interest in studying the Cholesky factorization was the solution of linear systems of equations Let’s first look at things through the lense of the recursive form of the Cholesky factorization. The first step of the factorization was to form the block Cholesky factorization<\/p>\n\n\n\n <\/span> <\/span> Suppose that we choose the top left To reduce computational cost, we ask the provocative question: What if we simply didn’t factorize the Schur complement? Observe that we can write the block Cholesky factorization as a sum of two terms<\/p>\n\n\n\n (2) <\/span> <\/span> We can use the first term of this sum as a rank- <\/span> <\/span> is the column Nystr\u00f6m approximation<\/em> (Nys) to <\/span> <\/span> Observe that the approximation Summing things up:<\/p>\n\n\n\n If we perform a partial<\/em> Cholesky factorization on a positive (semi)definite matrix, we obtain a low-rank approximation known as the column Nystr\u00f6m approximation. The residual of this approximation is the Schur complement, padded by rows and columns of zeros.<\/p><\/blockquote>\n\n\n\n The idea of obtaining a low-rank approximation from a partial matrix factorization is very common in matrix computations. Indeed, the optimal low-rank approximation to a real symmetric matrix is obtained by truncating its eigenvalue decomposition and the optimal low-rank approximation to a general matrix is obtained by truncating its singular value decomposition<\/a>. While the column Nystr\u00f6m approximation is not the optimal rank- The column Nystr\u00f6m approximation is formed from only Unfortunately there’s not a free lunch here. The column Nystr\u00f6m is only a good low-rank approximation if the Schur complement has small entries. In general, this need not be the case. Fortunately, we can improve our situation by means of pivoting<\/a><\/em>.<\/p>\n\n\n\n Our iterative Cholesky algorithm first performs elimination using the entry in position Obtaining good<\/em> column Nystr\u00f6m approximations requires identifying the a good choice for the pivots<\/em> to reduce the size of the entries of the Schur complement at each step of the algorithm. With general pivot selection, an iterative algorithm for computing a column Nystr\u00f6m approximation by partial Cholesky factorization proceeds as follows: Initialize an This procedure results in the Nystr\u00f6m approximation (Nys) with column set <\/span> <\/span> The pivoted Cholesky steps 1\u20133 requires updating the entire matrix How should we choose the pivots? Two simple strategies immediately suggest themselves:<\/p>\n\n\n\n <\/span> <\/span> The greedy strategy often (but not always) works well, and the uniformly random approach can work surprisingly well<\/a> if the matrix Despite often working fairly well, both the uniform and greedy schemes can fail significantly, producing very low-quality approximations. My research (joint with Yifan Chen<\/a>, Joel A. Tropp<\/a>, and Robert J. Webber<\/a>) has investigated a third strategy<\/a> striking a balance between these two approaches:<\/p>\n\n\n\n <\/span> <\/span> Our paper<\/a> provides theoretical analysis and empirical evidence showing that this diagonally-weighted random pivot selection (which we call randomly pivoted Cholesky<\/em> aka RPCholesky) performs well at approximating all positive semidefinite matrices In this post, we’ve seen that a column Nystr\u00f6m approximation can be obtained from a partial Cholesky decomposition. The residual of the approximation is the Schur complement. I hope you agree that this is a very nice connection between these three ideas. But beyond its mathematical beauty, why do we care about this connection? Here are my reasons for caring:<\/p>\n\n\n\n On a broader level, our tale of Nystr\u00f6m, Cholesky, and Schur demonstrates that there are rich connections between low-rank approximation and (partial versions of) classical matrix factorizations like LU<\/a> (with partial pivoting<\/a>), Cholesky<\/a>, QR<\/a>, eigendecomposition<\/a>, and SVD<\/a> for to full-rank matrices. These connections can be vital to analyzing low-rank approximation algorithms and developing improvements.<\/p>\n","protected":false},"excerpt":{"rendered":" In the last post, we looked at the Nystr\u00f6m approximation to an positive semidefinite (psd) matrix . A special case was the column Nystr\u00f6m approximation, defined to be (Nys) where identifies a subset of columns of . Provided , this allowed us to approximate all entries of the matrix using only the entries inRead more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[8],"tags":[],"class_list":["post-1370","post","type-post","status-publish","format-standard","hentry","category-low-rank-approximation-toolbox"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1370","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/comments?post=1370"}],"version-history":[{"count":18,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1370\/revisions"}],"predecessor-version":[{"id":1391,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1370\/revisions\/1391"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=1370"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=1370"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=1370"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}![\"\[\hat{A} = A(:,S) \, A(S,S)^{-1} \, A(S,:), \]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-cf768b8c5fedfd543a62070abae10bfe_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n![\"S = \{s_1,\ldots,s_k\} \subseteq \{1,2,\ldots,N\}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-e4e9d34e5b6acc4ed7c15b622e8fced3_l3.png\" "\\"Rendered")
identifies a subset of ![\"k\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-f65286b751f121928913d4aa91d94ee9_l3.png\" "\\"Rendered")
columns of . Provided ![\"k\ll N\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-c3709ff2e214db7f778a5f62fdbb9a42_l3.png\" "\\"Rendered")
, this allowed us to approximate all ![\"N^2\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d29e4554584afc4c0a6a384cedf37df8_l3.png\" "\\"Rendered")
entries of the matrix using only the ![\"kN\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-930a2318d303a9c8f951fe96800d8bd1_l3.png\" "\\"Rendered")
entries in columns ![\"s_1,\ldots,s_k\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-dce009c30c096fdcc766edae34514199_l3.png\" "\\"Rendered")
of , a huge savings of computational effort!<\/p>\n\n\n\n?<\/li>![\"A - \hat{A}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-a16c5891de2142e585fc70537203624d_l3.png\" "\\"Rendered")
of the approximation?<\/li><\/ul>\n\n\n\nCholesky: Solving Linear Systems<\/h2>\n\n\n\n
![\"Ax = b\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-bb1076898a88ac5d6a1c74d7932dd5fc_l3.png\" "\\"Rendered")
, where is a real invertible matrix and ![\"b\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-3ab6f6c64ce17f786a81f8cc8fdcec90_l3.png\" "\\"Rendered")
is a vector of length ![\"N\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-ec1be7928bcd319c17c1931fbb8059ac_l3.png\" "\\"Rendered")
. The standard way of doing this in modern practice (at least for non-huge matrices <\/a>) is by means of Gaussian elimination\/LU factorization<\/a>. We factor the matrix as a product ![\"A = LU\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-90c54eed56be2ab59ccb9688df92a1b7_l3.png\" "\\"Rendered")
of a lower triangular matrix ![\"L\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d99fc63db67b7ceb9f44b2bc4b03bc89_l3.png\" "\\"Rendered")
and an upper triangular matrix ![\"U\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-c7480669f4fca8a251671d27137d0b09_l3.png\" "\\"Rendered")
.2<\/a><\/sup>To make this accurate, we usually have to reorder the rows<\/a> of the matrix as well. Thus, we actually compute a factorization ![\"PA = LU\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-717b63dad015b09681392d4780ef8b0d_l3.png\" "\\"Rendered")
where ![\"P\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-e0b7f55ebbc9bb715e8b20ffdc459df7_l3.png\" "\\"Rendered")
is a permutation matrix and and are triangular.<\/span> The system is solved by first solving ![\"Ly = b\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d0cceead4f14e3035cafe7680ff4ce67_l3.png\" "\\"Rendered")
for ![\"y\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-7506eeeff09aad3bcf6b7259302df451_l3.png\" "\\"Rendered")
and then ![\"Ux = y\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-a5e82fd946fce03e6fade382ac3ec3dd_l3.png\" "\\"Rendered")
for ![\"x\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-b312d649591164b7149ed0756f694a76_l3.png\" "\\"Rendered")
; the triangularity of and make solving the associated systems of linear equations easy<\/a>.<\/p>\n\n\n\n, a simplification is possible. In this case, one can compute an LU factorization where the matrices and are transposes of each other, ![\"U = L^\top\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d5f18a8d1e3228d9ef5a17c5bca84aeb_l3.png\" "\\"Rendered")
. This factorization ![\"A = LL^\top\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-1e3f8839bf702fcaa4e7107093fa64a4_l3.png\" "\\"Rendered")
is known as a Cholesky factorization<\/a><\/strong> of the matrix .<\/p>\n\n\n\n to be complex-valued. For a complex-valued positive definite matrix<\/a> , its Cholesky decomposition takes the form ![\"A = LL^*\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-77bcb20c73da17bd82eb027d35d449ef_l3.png\" "\\"Rendered")
, where is again a lower triangular matrix. All that has changed is that the transpose ![\"{}^\top\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-0891b36c744e320081b679d0408bd179_l3.png\" "\\"Rendered")
has been replaced by the conjugate transpose<\/em><\/a> ![\"{}^*\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-c15a9e678f5eb8ef071c71b328820967_l3.png\" "\\"Rendered")
. We shall work with the more general complex case going forward, though the reader is free to imagine all matrices as real and interpret the operation as ordinary transposition if they so choose.<\/p>\n\n\n\nSchur: Computing the Cholesky Factorization<\/h2>\n\n\n\n
in block form<\/a> as<\/p>\n\n\n\n![\"\[A = \twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{22}}. \]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-24b23cd2ac7afc3e7c6acfb7fd8f8220_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n, factoring as a product of matrices which are only block triangular<\/em>. Then, we’ll “upgrade” this block factorization into a full Cholesky factorization.<\/p>\n\n\n\n by ![\"A_{21}A_{11}^{-1}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-aeb2a3885c54ebc2ecbe8d7fb870cd1a_l3.png\" "\\"Rendered")
and subtracting this from the second block row introduces a matrix of zeros into the bottom left block of . (The matrix ![\"A_{11}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-92456eb021fa3f381ff62422c0c24281_l3.png\" "\\"Rendered")
is a principal submatrix<\/em><\/a> of and is therefore guaranteed to be positive definite and thus invertible.3<\/a><\/sup>To directly see is positive definite, for instance, observe that since is positive definite, ![\"x^* A_{11}x = \twobyone{x}{0}^* A\twobyone{x}{0} > 0\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-2db932e5b62092a4bf58ee2122e892e0_l3.png\" "\\"Rendered")
for every nonzero vector .<\/span>) In matrix language,<\/p>\n\n\n\n![\"\[\twobytwo{I}{0}{-A_{21}A_{11}^{-1}}{I}\twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{22}} = \twobytwo{A_{11}}{A_{12}}{0}{A_{22} - A_{21}A_{11}^{-1}A_{12}}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-4b8f6deed9a1dd928c76d311ab88fd54_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n on the left-hand side of this equation by multiplying by<\/p>\n\n\n\n![\"\[\twobytwo{I}{0}{-A_{21}A_{11}^{-1}}{I}^{-1} = \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I}\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-be1610a1241e5c5bfd2c7dc04653e712_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n![\"\[A = \twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{22}} = \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I} \twobytwo{A_{11}}{A_{12}}{0}{A_{22} - A_{21}A_{11}^{-1}A_{12}}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-e004e13e06bf86fd8707805ce21f9071_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n into block triangular pieces, but these pieces are not (conjugate) transposes of each other. Thus, to make this equation more symmetrical, we can further decompose<\/p>\n\n\n\n![\"\[A = \twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{22}} = \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I} \twobytwo{A_{11}}{0}{0}{A_{22} - A_{21}A_{11}^{-1}A_{12}} \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I}^*. \]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-7e3a03d3b1e4421e71d17e2abba7ef97_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n taking the form ![\"A = LDL^*\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-0113ac536b3a0e5fd8dc447245849fa1_l3.png\" "\\"Rendered")
, where is a block lower triangular matrix and ![\"D\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-ab69a4a7716bbf890e5f604a06fd1f13_l3.png\" "\\"Rendered")
is a block diagonal matrix.<\/p>\n\n\n\n![\"\[S = A_{22} - A_{21}A_{11}^{-1}A_{12}. \]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-55e08632ed61d20591817565b66a7e09_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n![\"A_{22} - A_{21}A_{11}^{-1}A_{12}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-4ebb04dcb25be8e958d054cb126b53a0_l3.png\" "\\"Rendered")
in matrix theory makes it deserving of its special name, the Schur complement. To us for now, the Schur complement is just the matrix appearing in the bottom right corner of our block Cholesky factorization.<\/p>\n\n\n\n![\"A=\twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{22}}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-899e76844b6850c232b5bc8b080b66ef_l3.png\" "\\"Rendered")
is positive (semi)definite, then the Schur complement ![\"S = A_{22} - A_{21}A_{11}^{-1}A_{12}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-0ec78d555a1708661f29d2383c47c925_l3.png\" "\\"Rendered")
is positive (semi)definite.<\/p><\/blockquote>\n\n\n\n and ![\"S\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-b89471a11dd7fd85184a38ddb3ea9145_l3.png\" "\\"Rendered")
are positive definite.<\/p>\n\n\n\n![\"1\times 1\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-4ac03eedd19f9d68a98d70f3ba365091_l3.png\" "\\"Rendered")
matrix is ![\"A = A^{1/2} \cdot A^{1/2}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-9849ea39682b82f9ba5877fc5309e80d_l3.png\" "\\"Rendered")
.<\/span> computing Cholesky factorizations of the diagonal blocks<\/p>\n\n\n\n![\"\[A_{11} = L_{11}^{\vphantom{*}}L_{11}^*, \quad S = L_{22}^{\vphantom{*}}L_{22}^*.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-fcde81f118ea8516883b47d07a35b2b0_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n![\"LDL^*\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-a5ba506666ff7645c5b25148fae3a8a8_l3.png\" "\\"Rendered")
factorization (1) and simplifying gives a Cholesky factorization, as desired:<\/p>\n\n\n\n![\"\[A = \twobytwo{L_{11}}{0}{A_{21}^{\vphantom{*}}(L_{11}^{*})^{-1}}{L_{22}}\twobytwo{L_{11}}{0}{A_{21}^{\vphantom{*}}(L_{11}^{*})^{-1}}{L_{22}}^* =: LL^*.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-e22f8c5665f2154819e960cd6c6b32d3_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n!<\/p>\n\n\n\n to be of size ), our recursive Cholesky algorithm becomes the following iterative algorithm: Initialize to be the zero matrix. For ![\"j = 1,2,3,\ldots,N\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-60b5dafd0752ba10c33f4406f25b3006_l3.png\" "\\"Rendered")
, perform the following steps:<\/p>\n\n\n\n.<\/strong> Set the ![\"j\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-fe087f8cefab0bcb3270609914ada26c_l3.png\" "\\"Rendered")
th column of : ![\"\[L(j:N,j) \leftarrow A(j:N,j)/\sqrt{a_{jj}}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-9d8a1752f0ae4613f71767ce87b9560a_l3.png\" "\\"Rendered")
<\/p><\/li>.<\/strong> Update the bottom right portion of to be the Schur complement: ![\"\[A(j+1:N,j+1:N)\leftarrow A(j+1:N,j+1:N) - \frac{A(j+1:N,j)A(j,j+1:N)}{a_{jj}}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d8090111efb4d1beb891cac29c04a598_l3.png\" "\\"Rendered")
<\/p><\/li><\/ol>\n\n\n\nNystr\u00f6m: Using Cholesky Factorization for Low-Rank Approximation<\/h2>\n\n\n\n
for a positive definite matrix . We can also use the Cholesky factorization for a very different task, low-rank approximation.<\/p>\n\n\n\n![\"\[A = \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I} \twobytwo{A_{11}}{0}{0}{A_{22} - A_{21}A_{11}^{-1}A_{12}} \twobytwo{I}{0}{A_{21}A_{11}^{-1}}{I}^*.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-7d369f0415a2b62af274f183e8828853_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n block to be of size ![\"k\times k\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-979373ae59303770cd68d74797bc73f9_l3.png\" "\\"Rendered")
, where is much smaller than . The most expensive part of the Cholesky factorization will be the recursive factorization of the Schur complement ![\"A_{22} - A_{21}A_{11}^{-1} A_{12}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-9676b8dbd72a8d1166734af146a1fe48_l3.png\" "\\"Rendered")
, which is a large matrix of size ![\"(N-k)\times (N-k)\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-36c509a28ec6939da11d987fe90080a2_l3.png\" "\\"Rendered")
.<\/p>\n\n\n\n![\"\[A = \twobyone{I}{A_{21}{A_{11}^{-1}}} A_{11} \twobyone{I}{A_{22}{A_{11}^{-1}}}^* + \twobytwo{0}{0}{0}{A_{22}-A_{21}A_{11}^{-1}A_{12}}. \]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-eb7264120d2157941336da3dcbbcc1d8_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n approximation to the matrix <\/strong>. The low-rank approximation, which we can write out more conveniently as<\/p>\n\n\n\n![\"\[\hat{A} = \twobyone{A_{11}}{A_{21}} A_{11}^{-1} \onebytwo{A_{11}}{A_{12}} = \twobytwo{A_{11}}{A_{12}}{A_{21}}{A_{21}A_{11}^{-1}A_{12}},\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-8ec34b04edefd156d8de3711608998b2_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n associated with the index set ![\"S = \{1,2,3,\ldots,k\}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-f45c109a9719ae1d7d5618685b3d9eb8_l3.png\" "\\"Rendered")
and is the final of our three titular characters. The residual of the Nystr\u00f6m approximation is the second term in (2), which is none other than the Schur complement (Sch), padded by rows and columns of zeros:<\/p>\n\n\n\n![\"\[A - \hat{A} = \twobytwo{0}{0}{0}{A_{22}-A_{21}A_{11}^{-1}A_{12}}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-f8dccb2f8fd446335886ec9bdabc85da_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n![\"\hat{A}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-a9463bba1bf140207248e89a02db1929_l3.png\" "\\"Rendered")
is obtained from the process of terminating a Cholesky factorization midway through algorithm execution, so we say that the Nystr\u00f6m approximation results from a partial<\/em> Cholesky factorization of the matrix .<\/p>\n\n\n\n approximation to (though it does satisfy a weaker notion of optimality, as discussed in this previous post<\/a>), it has a killer feature not possessed by the optimal approximation:<\/p>\n\n\n\n columns<\/em> from the matrix . A column Nystr\u00f6m approximation approximates an matrix by only reading a fraction of its entries!<\/p><\/blockquote>\n\n\n\n![\"(1,1)\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-d324ea61edce693bff957ca1e663880c_l3.png\" "\\"Rendered")
followed by position ![\"(2,2)\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-8734ab1d9d809d2a2c9ae76a16ca761d_l3.png\" "\\"Rendered")
and so on. There’s no need to insist on this exact ordering of elimination steps. Indeed, at each step of the Cholesky algorithm, we can choose whichever diagonal position ![\"(j,j)\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-12d14fccb5184615fa205b44a0fce0c4_l3.png\" "\\"Rendered")
that we want to perform elimination. The entry we choose to perform elimination with is called a pivot<\/a><\/em>.<\/p>\n\n\n\n![\"N\times k\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-ba393b1539bb1c54c96343d6635ef3a6_l3.png\" "\\"Rendered")
matrix ![\"F\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-7bf5d1207baa8be58658ce9d3cf12043_l3.png\" "\\"Rendered")
to store the column Nystr\u00f6m approximation ![\"\hat{A} = FF^*\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-52eb639f2f3627d1a0801fcb6a7d1592_l3.png\" "\\"Rendered")
, in factored form. For ![\"j = 1,2,\ldots,k\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-e2799ec717f32c88edd3b14a4234552a_l3.png\" "\\"Rendered")
, perform the following steps:<\/p>\n\n\n\n![\"s_j\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-8c135c2dc6d5f751cde531932f5a063d_l3.png\" "\\"Rendered")
.<\/li>![\"F(:,j) \leftarrow A(:,s_j) / \sqrt{a_{s_js_j}}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-b1fdb907698e9174564ac35b7abc86e7_l3.png\" "\\"Rendered")
.<\/li>![\"A \leftarrow A - \frac{A(:,s_j)A(s_j,:)}{a_{s_js_j}}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-23132cecfc54e8fcbdf1adae27266c2e_l3.png\" "\\"Rendered")
.<\/li><\/ol>\n\n\n\n![\"S = \{s_1,\ldots,s_k\}\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-0156fa337c22e7e307bf93c47439dba9_l3.png\" "\\"Rendered")
:<\/p>\n\n\n\n![\"\[\hat{A} = FF^* = A(:,S) \, A(S,S)^{-1} \, A(S,:).\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-f9ac31fa249d889cd909be1d340abd58_l3.png\" "\\"Rendered")
<\/p><\/p>\n\n\n\n at every step. With a little more cleverness, we can optimize this procedure to only update the entries of the matrix we need to form the approximation . See Algorithm 2 in this preprint of my coauthors and I<\/a> for details.<\/p>\n\n\n\n, select uniformly at random from among the un-selected pivot indices.<\/li>, select according to the largest diagonal entry of the current residual : ![\"\[s_j = \argmax_{1\le k\le N} a_{kk}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-42a58dd2bd6824be5efdc930fc5356df_l3.png\" "\\"Rendered")
<\/p><\/li><\/ul>\n\n\n\n is “incoherent”, with all rows and columns of the matrix possessing “similar importance”.<\/p>\n\n\n\n, select at random according to the probability weights based on the current diagonal of the matrix ![\"\[\mathbb{P} \{ s_j = k \} = \frac{a_{kk}}{\operatorname{tr} A}.\]\"](\"https:\/\/www.ethanepperly.com\/wp-content\/ql-cache\/quicklatex.com-36ca3921c45ffb034e0cfa8adf85b982_l3.png\" "\\"Rendered")
<\/p><\/li><\/ul>\n\n\n\n, even those for which uniform random and greedy pivot selection fail. The success of this approach can be seen in the following examples (from Figure 1 in the paper<\/a>), which shows RPCholesky can produce much smaller errors than the greedy and uniform methods.<\/p>\n\n\n\n![\"\"](\"https:\/\/www.ethanepperly.com\/wp-content\/uploads\/2022\/10\/smile_trace-1024x768.png\")
<\/figure>\n\n\n\n![\"\"](\"https:\/\/www.ethanepperly.com\/wp-content\/uploads\/2022\/10\/outliers_trace-1024x768.png\")
<\/figure>\n<\/figure>\n\n\n\nConclusions<\/h2>\n\n\n\n