I've implemented the version of Aasen's algorithm described in the book Matrix Computations 4th Edition. The version there doesn't have pivoting. The book's description of how to add pivoting is a bit sketchy for me. Can someone illustrate with actual code how to add pivoting?

I should add that my implementation has support for Hermitian matrices, unlike the one in the book.

from sympy import *

def aasen(A):
    """Aasen without pivoting, but with support for Hermitian matrices."""
    n = A.rows
    alpha = zeros(n,1)
    beta = zeros(n,1)
    L = eye(n)
    h = zeros(n,1)
    v = zeros(n,1)
    for j in range(n):
        if j == 0:
            alpha[0] = A[0,0]
            v[1:n,0] = A[1:n,0]
            h[0] = beta[0].conjugate() * L[j,1].conjugate()
            for k in range(1,j):
                h[k] = beta[k-1] * L[j, k-1].conjugate() + alpha[k] * L[j,k].conjugate() + beta[k].conjugate() * L[j,k+1].conjugate()
            h[j] = A[j,j] - L[j,0:j].dot(h[0:j])
            alpha[j] = h[j] - beta[j-1] * L[j,j-1].conjugate()
            v[j+1:n,0] = A[j+1:n,j] - L[j+1:n,0:j+1] * Matrix(h[0:j+1])
        if j <= n-2:
            beta[j] = v[j+1]
        if j <= n - 3:
            L[j+2:n,j+1] = Matrix(v[j+2:n]) / v[j+1]
    T = zeros(n,n)
    for i in range(n):
        for j in range(n):
            if i == j+1:
                T[i,j] = beta[j]
            elif i == j:
                T[i,j] = alpha[j]
            elif j == i+1:
                T[i,j] = beta[i].conjugate()
    return L, T
  • 1
    $\begingroup$ Hello @ogogmad! It seems that after my answer on MO suggesting these algorithms you went all-in on implementing Aasen's. From what I understand, it is not much in use today: Bunch-Kaufman and Bunch-Parlett LDL are much more common; that is what Lapack's dsytrf uses, for instance. (And, in turn, every other software package basically just wraps Lapack). So I think you'd have a much easier time finding information on these LDL variants instead. $\endgroup$ Dec 5 '20 at 13:50
  • $\begingroup$ @FedericoPoloni I'm looking to produce or find Python code for any decomposition algorithm for Hermitian matrices. It doesn't really matter which one it is. So far, I'm nearly there when it comes to Aasen $\endgroup$
    – ogogmad
    Dec 5 '20 at 14:10
  • 1
    $\begingroup$ Though I highly respect your initiative in implement Aasen's method independently, I agree with Federico that you will have an easier time with this if you imitated the algorithms/vetted implementations already present in LAPACK. (Even if you're working beyond its s/d/c/z number types, the code would be basically the same). It would also have the advantage that if/when you do need to work with "plain" numbers, you could have your code delegate to numpy/LAPACK without much change in API. $\endgroup$ Dec 5 '20 at 18:32
  • 1
    $\begingroup$ Also, the complex Hermitian LDL^H decomposition is covered by zhetf2 / zhetrf / zhetrs family in LAPACK. (The zsytf2 / zsytrf / zsytrs family is for complex symmetric LDL^T). $\endgroup$ Dec 5 '20 at 18:35
  • 2
    $\begingroup$ @FedericoPoloni I'll try either Bunch-Parlett and/or Bunch-Kaufman $\endgroup$
    – ogogmad
    Dec 6 '20 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.