# Matrix diagonalization (linear algebra)

**URL:** <https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405>\
**Category:** New to Julia\
**Created:** [January 3, 2025, 9:36am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405 "2025-01-03T09:36:07Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 3, 2025, 9:36am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/1 "2025-01-03T09:36:07Z")

</div>

Hello everyone, I have a problem related to using Linear algebra package to help me diagonalize a matrix. Normally, one would use function like `eig()`. But as I found out that the returned eigenvalues are sorted automatically, so I am wondering if there is way I can keep the original order of my diagonal elements, because they actually have subindex which is ordered in my original matrix.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 3, 2025, 9:45am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/2 "2025-01-03T09:45:05Z")

</div>

Does your initial matrix have a special structure? Triangular?

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 3, 2025, 9:53am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/3 "2025-01-03T09:53:27Z")

</div>

Actually, right now it is 100\*100 matrix, so it is not big enough to consider optimization for now. Just testing the basic function.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 3, 2025, 9:59am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/4 "2025-01-03T09:59:46Z")

</div>

I was asking because I don’t understand why do you expect the order of your diagonal coefficients to be respected. The diagonal coefficients of an arbitrary matrix are not its eigenvalues, but you probably know that of course.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 3, 2025, 10:39am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/5 "2025-01-03T10:39:55Z")

</div>

It is not a special matrix. The reason is that each diagonal element is actually related a physical index, so if the output eigenvalue don’t respect that, then it will ruin the correspondence.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [January 3, 2025, 11:10am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/6 "2025-01-03T11:10:29Z")

</div>

Hm, that is not so easy, if you for example look at [Diagonalizable matrix - Wikipedia](https://en.wikipedia.org/wiki/Diagonalizable_matrix) – so assuming your matrix _is_ diagonalisable

P and D in A = P^{-1}DP are not unique in the following sense: If you permute the columns of P the elements in (on the diagonal of) D permute the same way.

So with that in mind there is nothing to respect, because there is no relation between the diagonal of A and the eigenvalues (=diagonal of D), you can chose _any_ permutation of the columns of P and obtain _any_ order.

`eig` had the convention to therefore order increasingly, which is still not unique for P, if eigenvalues repeat.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 3, 2025, 11:20am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/7 "2025-01-03T11:20:05Z")

</div>

Thanks, I can see the difficulty now, then I guess it means I have to use another approach to diagonalize the matrix.

---

<div class="post-metadata">

**Author:** ![Maucejo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maucejo/32/39090_2.png) [@Maucejo](https://discourse.julialang.org/u/Maucejo)\
**Post date:** [January 3, 2025, 1:35pm UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/8 "2025-01-03T13:35:43Z")

</div>

I suggest you to have a look at the documentation of the `eigen` function of `LinearAlgebra` and the `sortby` keyword.

> **[Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.eigen)**
>
> Documentation for The Julia Language.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 3, 2025, 1:58pm UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/9 "2025-01-03T13:58:52Z")

</div>

> [@fanzhy171](#):
>
> then I guess it means I have to use another approach to diagonalize the matrix.

There is no other approach. Diagonalizing a matrix is _synonymous_ with finding eigenvalues and eigenvectors: you are finding a change of basis (eigenvectors) that makes the matrix diagonal (the eigenvalues). I think you need to re-think your whole conception of what diagonalization means and how it relates to your problem.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 2:28am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/10 "2025-01-04T02:28:53Z")

</div>

Thanks, I read the documents before, but I don’t see any example on how to design a custom `sortby` function. Also As @kellertuer was mentioned, in the numerical procedure I am afraid there is no way that I can keep track of every diagonal elements.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 2:38am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/11 "2025-01-04T02:38:41Z")

</div>

Actually I am thinking use a way which is called _flow equation_ method, which will make diagonalization process becomes solving coupled differential equations. The amazing part of this method is that you can keep track of the flow of the diagonal and off-diagonal part of your matrix, at the end of the day you can see all the off-diagonal coefficients are gone _guaranteed_ (Btw, this method is implemented to solve one-body quantum problems, the Hamiltonian is our matrix). The reference is:

[https://onlinelibrary.wiley.com/doi/10.1002/andp.19945060203](https://onlinelibrary.wiley.com/doi/10.1002/andp.19945060203)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 4, 2025, 3:26am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/12 "2025-01-04T03:26:22Z")

</div>

> [@fanzhy171](#):
>
> Actually I am thinking use a way which is called _flow equation_ method, which will make diagonalization process becomes solving coupled differential equations.

Of course there are _many_ algorithms to find a matrix diagonalization — but at the end of the day they are **all computing equivalent things** : the diagonal matrix _must_ consist of the eigenvalues (in some order), and the change of basis to make it diagonal must consist of eigenvectors.

The ordering of the eigenvalues can always be permuted arbitrarily by permuting the eigenvectors correspondingly. Because of this, it’s not clear what connection you are hoping to maintain between the diagonal entries of your original matrix (?) and the ordering of the eigenvalues (the diagonal entries of the diagonalized matrix).

Maybe the flow method yields the ordering you want? But if there is some symmetry that you are trying to preserve I would try to express that more explicitly algebraically, to give you more freedom in exploiting state-of-the-art eigensolver algorithms. e.g. if you have a group of symmetry operators that commutes with your Hamiltonian then you can use it to categorize and partially order the eigenvectors (according to the irreps of the symmetry group).

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 4:41am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/13 "2025-01-04T04:41:25Z")

</div>

To make it clear, my original matrix has momentum indices, and the diagonal elements represents a dispersion relationship, there are off-diagonal couplings in this case. I want to see after this diagonalization, what will my band dispersion get renormalized, so if the eigenvalue solver messed up my coherence, meaning the eigenvalue’s indices are not corresponding to my original indices, then I can not figuring that out easily. And I think Flow equation method intrinsically protect this property. Hope this can make my problem more clear, and thx for your question and response.

---

<div class="post-metadata">

**Author:** ![linwaytin](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@linwaytin](https://discourse.julialang.org/u/linwaytin)\
**Post date:** [January 4, 2025, 4:58am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/14 "2025-01-04T04:58:48Z")

</div>

You said “my original matrix has momentum indices”. Do you mean you discretize a Hamiltonian H(k,k') to get a matrix H\_{i,j} with i,j corresponding to some momenta?

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 5:00am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/15 "2025-01-04T05:00:19Z")

</div>

Yes, all momentum indices not real space.

---

<div class="post-metadata">

**Author:** ![linwaytin](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@linwaytin](https://discourse.julialang.org/u/linwaytin)\
**Post date:** [January 4, 2025, 5:37am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/16 "2025-01-04T05:37:55Z")

</div>

Then the eigenvalue problem is H\_{k,k'}u\_{k'}=Eu\_{k'}. As many people have said in the above, the eigenvalue E in general has no relation to the original k. You need to introduce additional structure to “keep” the original indices.

I think the confusion is from perturbation theory, where you can discuss the perturbed eigenvalues and thus you can keep the original index. Maybe this is related to the paper you mentioned.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 6:52am UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/17 "2025-01-04T06:52:09Z")

</div>

Yes, this can be understood as a perturbation. But the method I mentioned before doesn’t necessarily need perturbation assumptions.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 4, 2025, 12:15pm UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/18 "2025-01-04T12:15:17Z")

</div>

I think what you might be able to do is using the eigenvectors to relate the eigenvalues to their “original” position by considering the overlaps between the original basis vectors and the eigenvectors. Since the original basis vectors are just vectors where one index is 1 and the rest are zeros, you can simply take the absolute values of the components of each eigenvector and the index of the largest value gives you the “original” basis index (which is the vector with the greatest overlap). If the “perturbations” are reasonably small, this should work well. However if there are degeneracies in your original matrix, then the corresponding states could hybridize and then your question is not quite well-defined I think.  
If you are using a flow equation method, it might make sense to keep track of these correspondences along the flow. IIUC, the eigenvalues/eigenvectors are deformed continuously along the flow (well in the absence of degeneracies at least) so this approach should work well.

---

<div class="post-metadata">

**Author:** ![fanzhy171](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fanzhy171/32/213949_2.png) [@fanzhy171](https://discourse.julialang.org/u/fanzhy171)\
**Post date:** [January 4, 2025, 1:40pm UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/19 "2025-01-04T13:40:54Z")

</div>

Thanks for your response and suggestion. In my problem, the effect of the off-diagonal coupling should generally not work in a perturbative sense. A easy version of my problem is to consider a plane wave scattering off from a negative potential, and bound state solution can’t be generated from the perturbation type calculation. Yes you are absolutely right on the idea of flow equation method.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 4, 2025, 4:35pm UTC](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405/20 "2025-01-04T16:35:14Z")

</div>

I wasn’t able to look at the paper you cited to see exactly what the flow equations are. That was partly because I lack institutional access to the online article, but also I doubt I have enough background in physics to find the explanations clear. Explanations in terms of physics, instead of just math, will do more to confuse me than to enlighten me.

With very limited understanding of what you are trying to do, what I have taken away from the discussion is that you want to diagonalize and end up with the eigenvalues and eigenvectors in the same order you would have gotten with the flow equations, without actually using the flow equations for the computation. That sounds hard to do unless you know something in advance about the ordering of eigenvalues your flow will give you.

With that in mind, in case you don’t know, I wanted to point out that there is prior work in the literature on numerical analysis on diagonalizing with flows. You might want to look at the paper [Deift, Nanda, and Tomei, “Ordinary Differential Equations and the Symmetric Eigenvalue Problem”](https://www.jstor.org/stable/2157167). Some flows do produce sorted eigenvalues in the sense that the only stable equilibrium point is a diagonal matrix with the eigenvalues sorted in descending order. Maybe there is some ordering imposed by your flow that you can exploit.

[Next page](https://discourse.julialang.org/t/matrix-diagonalization-linear-algebra/124405.md?page=2)