Bunch-Kaufman with rook pivoting

General Usage
1

I asked this question on Slack but didn’t get an answer. Could someone explain how to use the Bunch-Kaufman factorization with rook pivoting?

julia> A = [0.271242 0.657867 1.14586; 0.657867 1.20111 1.13588; 1.14586 1.13588 0.434028]
3×3 Array{Float64,2}:
 0.271242  0.657867  1.14586
 0.657867  1.20111   1.13588
 1.14586   1.13588   0.434028

julia> SA = Symmetric(A, :L);

julia> F = bunchkaufman(SA, true, check=true);  # rook pivoting = true

julia> F.L * F.D * F.L' ≈ A[F.p, F.p]  # this part appears to succeed with rook = false
false

julia> F.L * F.D * F.L' ≈ A
false

How are we to use the permutation to recover the original matrix? Thanks.

2

It looks like this is a bug:

and the rook-pivoted Bunch–Kaufman factorization is returning the wrong permutation.

Note, however, than most of the time you use the factorization object simply by x = F \ b to solve Ax=b, and the \ solver for rook-pivoted Bunch–Kaufman appears to be correct.

1 Like
3

I noticed that \ appears to return the correct value (presumably because it calls LAPACK directly instead of relying on the permutation). However, in my use case I specifically need access to L and D.

4

In case you haven’t already derived the fix from the LAPACK docs, I think the following may do the job:

# adapted from LinearAlgebra stdlib
function _ipiv2p_rook(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar) where T
    p = T[1:maxi;]
    uploL = uplo == 'L'
    i = uploL ? 1 : maxi
    # if uplo == 'U' we construct the permutation backwards
    @inbounds while 1 <= i <= length(v)
        vi = v[i]
        if vi > 0 # the 1x1 blocks
            p[i], p[vi] = p[vi], p[i]
            i += uploL ? 1 : -1
        else # the 2x2 blocks
            if uploL
                p[i], p[-vi] = p[-vi], p[i]
                vp = -v[i+1]
                p[i + 1], p[vp] = p[vp], p[i + 1]
                i += 2
            else # 'U'
                p[i], p[-vi] = p[-vi], p[i]
                vp = -v[i-1]
                p[i - 1], p[vp] = p[vp], p[i - 1]
                i -= 2
            end
        end
    end
    return p
end
function permvec(B::BunchKaufman{T}) where {T}
    B.rook || return B.p
    n = size(B,1)
    return _ipiv2p_rook(getfield(B, :ipiv), n, getfield(B, :uplo))
end

Then

F = bunchkaufman(A,rook)
p = permvec(F)
F.L * F.D * tfunc(F.L) ≈ A[p,p]

where tfunc is adjoint or transpose as appropriate.

5 Likes
5

A PR for Base to fix #32080 would be welcome.

1 Like
6

Many thanks @Ralph_Smith! That does the job and appears to match the LAPACK documentation.