{"id":1843,"date":"2024-06-25T19:28:53","date_gmt":"2024-06-25T19:28:53","guid":{"rendered":"https:\/\/www.ethanepperly.com\/?p=1843"},"modified":"2024-06-28T22:00:41","modified_gmt":"2024-06-28T22:00:41","slug":"neat-randomized-algorithms-randomized-cholesky-qr","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2024\/06\/25\/neat-randomized-algorithms-randomized-cholesky-qr\/","title":{"rendered":"Neat Randomized Algorithms: Randomized Cholesky QR"},"content":{"rendered":"\n

As a research area, randomized numerical linear algebra (RNLA) is as hot as ever. To celebrate the exciting work in this space, I’m starting a new series on my blog where I celebrate cool recent algorithms in the area. At some future point, I might talk about my own work in this series, but for now I’m hoping to use this series to highlight some of the awesome work being done by my colleagues.<\/em><\/p>\n\n\n\n

\n\n\n\n

Given a tall matrix \"A \in \real^{m\times n}\" with \"m\ge n\", its (economy-size<\/a>) QR factorization<\/a> is a decomposition of the form \"A = QR\", where \"Q \in \real^{m\times n}\" is a matrix with orthonormal columns<\/a> and \"R \in \real^{n\times n}\" is upper triangular<\/a>. QR factorizations are used to solve least-squares problems<\/a> and as a computational procedure to orthonormalize the columns of a matrix.<\/p>\n\n\n\n

Here’s an example in MATLAB, where we use QR factorization to orthonormalize the columns of a \"10^6\times 10^2\" test matrix. It takes about 2.5 seconds to run.<\/p>\n\n\n\n

>> A = randn(1e6, 1e2) * randn(1e2) * randn(1e2); % test matrix<\/mark>\n>> tic; [Q,R] = qr(A,\"econ\"); toc\nElapsed time is 2.647317 seconds.<\/code><\/pre>\n\n\n\n
The classical algorithm<\/a> for computing a QR factorization uses Householder reflectors<\/a> and is exceptionally numerically stable. Since \"Q\" has orthonormal columns, \"Q^\top Q = I\" is the identity matrix<\/a>. Indeed, this relation holds up to a tiny error for the \"Q\" computed by Householder QR:<\/p>\n\n\n\n
>> norm(Q'*Q - eye(1e2)) % || Q^T Q - I ||<\/mark>\nans =\n   7.0396e-14<\/code><\/pre>\n\n\n\n
The relative error \"\|A - QR\|/\|A\|\" is also small:<\/p>\n\n\n\n
>> norm(A - Q*R) \/ norm(A)\nans =\n   4.8981e-14<\/code><\/pre>\n\n\n\n
Here is an alternate procedure for computing a QR factorization, known as Cholesky QR:<\/p>\n\n\n\n
function [Q,R] = cholesky_qr(A)\n   R = chol(A'*A);\n   Q = A \/ R;       % Q = A * R^{-1}<\/mark>\nend<\/code><\/pre>\n\n\n\n
This algorithm works by forming \"A^\top A\", computing its (upper triangular) Cholesky decomposition<\/a> \"A^\top A = R^\top R\", and setting \"Q = AR^{-1}\". Cholesky QR is very fast, about \"5\times\" faster than Householder QR for this example:<\/p>\n\n\n\n
>> tic; [Q,R] = cholesky_qr(A); toc\nElapsed time is 0.536694 seconds.<\/code><\/pre>\n\n\n\n
Unfortunately, Cholesky QR is much less accurate and numerically stable than Householder QR. Here, for instance, is the value of \"\|Q^\top Q - I\|\", about ten million times larger than for Householder QR!:<\/p>\n\n\n\n
>> norm(Q'*Q - eye(1e2))\nans =\n   7.5929e-07<\/code><\/pre>\n\n\n\n
What’s going on? As we’ve discussed before<\/a> on this blog<\/a>, forming \"A^\top A\" is typically problematic in linear algebraic computations<\/a>. The “badness” of a matrix \"A\" is measured by its condition number<\/a><\/em>, defined to be the ratio of its largest and smallest singular values<\/a> \"\kappa(A) = \sigma_{\rm max}(A)/\sigma_{\rm min}(A)\". The condition number of \"A^\top A\" is the square of the condition number of \"A\", \"\kappa(A^\top A) = \kappa(A)^2\", which is at the root of Cholesky QR’s loss of accuracy. Thus, Cholesky QR is only appropriate for matrices that are well-conditioned<\/em>, having a small condition number \"\kappa(A) \approx 1\", say \"\kappa(A) \le 10\".<\/p>\n\n\n\n
The idea of randomized Cholesky QR<\/strong> is to use randomness to precondition<\/a><\/em> \"A\", producing a matrix \"R_1\" that \"B = AR_1^{-1}\" is well-conditioned. Then, since \"B\" is well-conditioned, we can apply ordinary Cholesky QR to it without issue. Here are the steps:<\/p>\n\n\n\n
    \n
  1. Draw a sketching matrix \"S\" of size \"2n\times m\"; see these<\/a> posts<\/a> of mine for an introduction to sketching.<\/li>\n\n\n\n
  2. Form the sketch \"SA\". This step compresses the very tall matrix \"m\times n\" to the much shorter matrix \"SA\" of size \"2n\times n\".<\/li>\n\n\n\n
  3. Compute a QR factorization \"SA = Q_1R_1\" using Householder QR. Since the matrix \"SA\" is small, this factorization will be quick to compute.<\/li>\n\n\n\n
  4. Form the preconditioned matrix \"B = AR_1^{-1}\".<\/li>\n\n\n\n
  5. Apply Cholesky QR to \"B\" to compute \"B = QR_2\".<\/li>\n\n\n\n
  6. Set \"R = R_2R_1\". Observe that \"A = BR_1 = QR_2R_1 = QR\", as desired.<\/li>\n<\/ol>\n\n\n\n
    MATLAB code for randomized Cholesky QR is provided below:1<\/a><\/sup>Code for the sparsesign subroutine can be found here<\/a>.<\/span>\n\n\n\n
    function [Q,R] = rand_cholesky_qr(A)\n   S = sparsesign(2*size(A,2),size(A,1),8); % sparse sign embedding<\/mark>\n   R1 = qr(S*A,\"econ\"); % sketch and (Householder) QR factorize<\/mark>\n   B = A \/ R1; % B = A * R_1^{-1}<\/mark>\n   [Q,R2] = cholesky_qr(B);\n   R = R2*R1;\nend<\/code><\/pre>\n\n\n\n
    Randomized Cholesky QR is still faster than ordinary Householder QR, about \"3\times\" faster in our experiment:<\/p>\n\n\n\n
    >> tic; [Q,R] = rand_cholesky_qr(A); toc\nElapsed time is 0.920787 seconds.<\/code><\/pre>\n\n\n\n
    Randomized Cholesky QR greatly improves on ordinary Cholesky QR in terms of accuracy and numerical stability. In fact, the size of \"\|Q^\top Q - I\|\" is even smaller than for Householder QR!<\/p>\n\n\n\n
    >> norm(Q'*Q - eye(1e2))\nans =\n   1.0926e-14<\/code><\/pre>\n\n\n\n
    The relative error \"\|A - QR\|/\|A\|\" is small, too! Even smaller than for Householder QR in fact:<\/p>\n\n\n\n
    >> norm(A - Q*R) \/ norm(A)\nans =\n   4.0007e-16<\/code><\/pre>\n\n\n\n
    Like many great ideas, randomized Cholesky QR was developed independently by a number of research groups. A version of this algorithm was first introduced in 2021 by Fan, Guo, and Lin<\/a>. Similar algorithms were investigated in 2022 and 2023 by Balabanov<\/a>, Higgins et al.<\/a>, and Melnichenko et al.<\/a> Check out Melnichenko et al.<\/a>‘s paper in particular, which shows very impressive results for using randomized Cholesky QR to compute column pivoted<\/em> QR factorizations<\/a>.<\/p>\n\n\n\n
    References: <\/strong>Primary references are A Novel Randomized XR-Based Preconditioned CholeskyQR Algorithm<\/a><\/em> by Fan, Guo, and Lin (2021); Randomized Cholesky QR factorizations<\/a><\/em> by Balabanov (2022); Analysis of Randomized Householder-Cholesky QR Factorization with Multisketching<\/a><\/em> by Higgins et al. (2023); CholeskyQR with Randomization and Pivoting for Tall Matrices (CQRRPT)<\/a><\/em> by Melnichenko et al. (2023). The idea of using sketching to precondition tall matrices originates in the paper A fast randomized algorithm for overdetermined linear least-squares regression<\/a><\/em> by Rokhlin and Tygert (2008).<\/p>\n","protected":false},"excerpt":{"rendered":"
    As a research area, randomized numerical linear algebra (RNLA) is as hot as ever. To celebrate the exciting work in this space, I’m starting a new series on my blog where I celebrate cool recent algorithms in the area. At some future point, I might talk about my own work in this series, but forRead 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":[6,15,12],"tags":[],"class_list":["post-1843","post","type-post","status-publish","format-standard","hentry","category-expository","category-neat-randomized-algorithms","category-sketching"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1843","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=1843"}],"version-history":[{"count":10,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1843\/revisions"}],"predecessor-version":[{"id":1867,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1843\/revisions\/1867"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=1843"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=1843"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=1843"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}