Fast `diag(A' * B * A)`

[Mason]

@Elrod

Also, FWIW, it may be faster to do

    D = Vector{promote_type(T,U)}(undef, size(A, 2))
    @tturbo for i ∈ eachindex(D)
        di = zero(eltype(D))
        for j ∈ axes(B, 1)
            for k ∈ axes(B, 2)
                di += conj(A[j, i]) * B[j, k] * A[k, i]
            end
        end
        D[i] = di
    end

Shaves off a couple hundred nanoseconds fore me, I’ll add it to the benchmark pile.

[e3c6]

If I add the assumption that B is Hermitian, would it help?

[mikmoore]

If B is Hermitian and positive definite then it can be factorized as B = C'*C. cholesky or sqrt are capable of producing such factorizations.

One can use the fact that diag(X'*X) == vec(sum(abs2, X; dims=1)). Setting X = C*A such that X'*X == A'*B*A can allow one to do this. But this still requires factorizing B = C'*C and then multiplying C*A. C is the same size as B. It may be triangular, which could save half the multiplications, but the cost of computing the factorization in the first place is likely larger than the savings. So ultimately this approach appears to be a net-loss in performance even under good circumstances.

One could also use a factorization B = C'*D*C (which does not require positive definiteness), but this has the same issues.

So the short answer is “no”, I don’t think B being Hermitian helps you to compute this any more quickly (barring exceptional structure). An exception would be if you already have B in some factorized form, i.e. that you computed it from C'*C for some C. Then you should just use the factors C directly in your calculations (which may spare you from ever needing to form B at all).