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

**URL:** <https://discourse.julialang.org/t/fast-diag-a-b-a/98216>\
**Category:** General Usage\
**Tags:** performance, linearalgebra\
**Created:** [May 2, 2023, 5:46pm UTC](https://discourse.julialang.org/t/fast-diag-a-b-a/98216 "2023-05-02T17:46:25Z")\
**Posts on this page:** 3\
**Page:** 2

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [May 2, 2023, 10:41pm UTC](https://discourse.julialang.org/t/fast-diag-a-b-a/98216/21 "2023-05-02T22:41:51Z")

</div>

> [@Elrod](#):
>
> Also, FWIW, it may be faster to do
> 
> ```julia
> 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.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 4, 2023, 7:14am UTC](https://discourse.julialang.org/t/fast-diag-a-b-a/98216/22 "2023-05-04T07:14:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 4, 2023, 3:32pm UTC](https://discourse.julialang.org/t/fast-diag-a-b-a/98216/23 "2023-05-04T15:32:22Z")

</div>

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).

[Previous page](https://discourse.julialang.org/t/fast-diag-a-b-a/98216.md?page=1)