{"id":1852,"date":"2024-07-02T16:04:48","date_gmt":"2024-07-02T16:04:48","guid":{"rendered":"https:\/\/www.ethanepperly.com\/?p=1852"},"modified":"2024-07-02T16:04:49","modified_gmt":"2024-07-02T16:04:49","slug":"neat-randomized-algorithms-randdiag-for-rapidly-diagonalizing-normal-matrices","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2024\/07\/02\/neat-randomized-algorithms-randdiag-for-rapidly-diagonalizing-normal-matrices\/","title":{"rendered":"Neat Randomized Algorithms: RandDiag for Rapidly Diagonalizing Normal Matrices"},"content":{"rendered":"\n

Consider two complex-valued square matrices \"A\in\complex^{n\times n}\" and \"B\in\complex^{n\times n}\". The first matrix \"A\" is Hermitian<\/a>, being equal \"A = A^*\" to its conjugate transpose<\/a> \"A^*\". The other matrix \"B\" is non-Hermitian, \"B \ne B^*\". Let’s see how long it takes to compute their eigenvalue decompositions in MATLAB:<\/p>\n\n\n\n

>> A = randn(1e3) + 1i*randn(1e3); A = (A+A')\/2;\n>> tic; [V_A,D_A] = eig(A); toc % Hermitian matrix<\/mark>\nElapsed time is 0.415145 seconds.\n>> B = randn(1e3) + 1i*randn(1e3);\n>> tic; [V_B,D_B] = eig(B); toc % non-Hermitian matrix<\/mark>\nElapsed time is 1.668246 seconds.<\/code><\/pre>\n\n\n\n
We see that it takes \"4\times\" longer to compute the eigenvalues of the non-Hermitian matrix \"B\" as compared to the Hermitian matrix \"A\". Moreover, the matrix \"V_A\" of eigenvectors for a Hermitian matrix \"A = V_AD_AV_A^{-1}\" is a unitary matrix<\/a>, \"V_A^*V_A = V_AV_A^* = I\".<\/p>\n\n\n\n
There are another class of matrices with nice eigenvalue decompositions, normal matrices<\/a>. A square, complex-valued matrix \"C\" is normal<\/a><\/em> if \"C^*C = CC^*\". The matrix \"V_C\" of eigenvectors for a normal matrix \"C = V_C D_C V_C^{-1}\" is also unitary, \"V_C^*V_C = V_CV_C^* = I\". An important class of normal matrices are unitary matrices themselves. A unitary matrix \"U\" is always normal since it satisfies \"U^*U = UU^* = I\".<\/p>\n\n\n\n
Let’s see how long it takes MATLAB to compute the eigenvalue decomposition of a unitary (and thus normal) matrix:<\/p>\n\n\n\n
>> U = V_A;                     % unitary, and thus normal, matrix<\/mark>\n>> tic; [V_U,D_U] = eig(U); toc % normal matrix<\/mark>\nElapsed time is 2.201017 seconds.<\/code><\/pre>\n\n\n\n
Even longer than it took to compute an eigenvalue decomposition of the non-normal matrix \"B\"! Can we make the normal eigenvalue decomposition closer to the speed of the Hermitian eigenvalue decomposition?<\/p>\n\n\n\n
Here is the start of an idea. Every square matrix \"C\" has a Cartesian decomposition<\/em>:
   <\/span>   <\/span>\"\[C = H + iS, \quad H = \frac{C+C^*}{2}, \quad S = \frac{C-C^*}{2i}.\]\"<\/p>We have written \"C\" as a combination of its Hermitian part<\/em> \"H\" and \"i\" times its skew-Hermitian part<\/em> \"S\". Both \"H\" and \"S\" are Hermitian matrices. The Cartesian decomposition of a square matrix is analogous to the decomposition of a complex number into its real and imaginary parts.<\/p>\n\n\n\n
For a normal matrix \"C\", the Hermitian and skew-Hermitian parts commute, \"HS = SH\". We know from matrix theory that commuting Hermitian matrices are simultaneously diagonalizable<\/a>, i.e., there exists \"Q\" such that \"H = QD_HQ^*\" and \"S = QD_SQ^*\" for diagonal matrices \"D_H\" and \"D_S\". Thus, given access to such \"Q\", \"C\" has eigenvalue decomposition 
   <\/span>   <\/span>\"\[C = Q(D_H+iD_S)Q^*.\]\"<\/p><\/p>\n\n\n\n
Here’s a first attempt to turn this insight into an algorithm. First, compute the Hermitian part \"H\" of \"C\", diagonalize \"H = QD_HQ^*\", and then see if \"Q\" diagonalizes \"C\". Let’s test this out on a \"2\times 2\" example:<\/p>\n\n\n\n
>> C = orth(randn(2) + 1i*randn(2)); % unitary matrix<\/mark>\n>> H = (C+C')\/2;                     % Hermitian part<\/mark>\n>> [Q,~] = eig(H);\n>> Q'*C*Q                            % check to see if diagonal<\/mark>\nans =\n  -0.9933 + 0.1152i  -0.0000 + 0.0000i\n   0.0000 + 0.0000i  -0.3175 - 0.9483i<\/code><\/pre>\n\n\n\n
Yay! We’ve succeeded at diagonalizing the matrix \"C\" using only a Hermitian eigenvalue decomposition. But we should be careful about declaring victory too early. Here’s a bad example:<\/p>\n\n\n\n
>> C = [1 1i;1i 1]; % normal matrix<\/mark>\n>> H = (C+C')\/2;\n>> [Q,~] = eig(H);\n>> Q'*C*Q           % oh no! not diagonal<\/mark>\nans =\n   1.0000 + 0.0000i   0.0000 + 1.0000i\n   0.0000 + 1.0000i   1.0000 + 0.0000i<\/code><\/pre>\n\n\n\n
What’s going on here? The issue is that the Hermitian part \"H = I\" for this matrix has a repeated eigenvalue. Thus, \"H\" has multiple different valid matrices of eigenvectors. (In this specific case, every<\/em> unitary matrix \"Q\" diagonalizes \"H\".) By looking at \"H\" alone, we don’t know which \"Q\" matrix to pick which also diagonalizes \"S\".<\/p>\n\n\n\n
He and Kressner developed a beautifully simple randomized<\/em> algorithm called RandDiag<\/a><\/em> to circumvent this failure mode. The idea is straightforward:<\/p>\n\n\n\n
    \n
  1. Form a random linear combination \"M = \gamma_1 H + \gamma_2 S\" of the Hermitian and skew-Hermitian parts of \"A\", with standard normal<\/a> random coefficients \"\gamma_1\" and \"\gamma_2\".<\/li>\n\n\n\n
  2. Compute \"Q\" that diagonalizes \"M\".<\/li>\n<\/ol>\n\n\n\n
    That’s it!<\/p>\n\n\n\n
    To get a sense of why He and Kressner’s algorithm works, suppose that \"H\" has some repeated eigenvalues and \"S\" has all distinct eigenvalues. Given this setup, it seems likely that a random linear combination of \"S\" and \"H\" will also have all distinct eigenvalues. (It would take a very special circumstances for a random linear combination to yield two eigenvalues that are exactly<\/em> the same!) Indeed, this intuition is a fact: With 100% probability,<\/a> \"Q\" diagonalizing a Gaussian random linear combination of simultaneously diagonalizable matrices \"H\" and \"S\" also diagonalizes \"H\" and \"S\" individually.<\/p>\n\n\n\n
    MATLAB code for RandDiag is as follows:<\/p>\n\n\n\n
    function Q = rand_diag(C)\n   H = (C+C')\/2; S = (C-C')\/2i;\n   M = randn*H + randn*S;\n   [Q,~] = eig(M);\nend<\/code><\/pre>\n\n\n\n
    When applied to our hard \"2\times 2\" example from before, RandDiag succeeds at giving a matrix that diagonalizes \"C\":<\/p>\n\n\n\n
    >> Q = rand_diag(C);\n>> Q'*C*Q\nans =\n   1.0000 - 1.0000i  -0.0000 + 0.0000i\n  -0.0000 - 0.0000i   1.0000 + 1.0000i<\/code><\/pre>\n\n\n\n
    For computing the matrix of eigenvectors for a \"1000\times 1000\" unitary matrix, RandDiag takes 0.4 seconds, just as fast as the Hermitian eigendecomposition did.<\/p>\n\n\n\n
    >> tic; V_U = rand_diag(U); toc\nElapsed time is 0.437309 seconds.<\/code><\/pre>\n\n\n\n
    He and Kressner’s algorithm is delightful. Ultimately, it uses randomness in only a small way. For most coefficients \"a_1,a_2 \in \real\", a matrix \"Q\" diagonalizing \"a_1 H + a_2 S\" will also diagonalize \"A = H+iS\". But, for any specific choice of \"a_1,a_2\", there is a possibility of failure. To avoid this possibility, we can just pick \"a_1\" and \"a_2\" at random. It’s really as simple as that.<\/p>\n\n\n\n
    References:<\/strong> RandDiag was proposed in A simple, randomized algorithm for diagonalizing normal matrices<\/a><\/em> by He and Kressner (2024), building on their earlier work in Randomized Joint Diagonalization of Symmetric Matrices<\/a><\/em> (2022) which considers the general case of using random linear combinations to (approximately) simultaneous diagonalize (nearly) commuting matrices. RandDiag is an example of a linear algebraic algorithm that uses randomness to put the input into “general position”; see Randomized matrix computations: Themes and variations<\/a><\/em> by Kireeva and Tropp (2024) for a discussion of this, and other, ways of using randomness to design matrix algorithms.<\/p>\n","protected":false},"excerpt":{"rendered":"
    Consider two complex-valued square matrices and . The first matrix is Hermitian, being equal to its conjugate transpose . The other matrix is non-Hermitian, . Let’s see how long it takes to compute their eigenvalue decompositions in MATLAB: We see that it takes longer to compute the eigenvalues of the non-Hermitian matrix as compared toRead 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],"tags":[],"class_list":["post-1852","post","type-post","status-publish","format-standard","hentry","category-expository","category-neat-randomized-algorithms"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1852","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=1852"}],"version-history":[{"count":7,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1852\/revisions"}],"predecessor-version":[{"id":1859,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1852\/revisions\/1859"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=1852"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=1852"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=1852"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}