# 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:** 1\
**Showing post:** 9

<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 2, 2023, 7:58pm UTC](https://discourse.julialang.org/t/fast-diag-a-b-a/98216/9 "2023-05-02T19:58:56Z")

</div>

> [@Mason](#):
>
> `D[i] = conj(A[j, i]) * B[j, k] * A[j, i]`

Sorry to not catch this all at once.

Shouldn’t this be `+=` and initialized to `D = zeros(V, size(A,2))`? Or maybe use a temporary variable `Di = zero(V)` to do the accumulation before saving to `D[i]`, if the compiler won’t hoist that on its own…

As written, the loop only does “work” on the final pair of `j,k`. The compiler may not be smart enough to have caught on, so it may not have a performance impact, but I’m suspicious that this function manages to be so much faster than even `B*A` despite that it should be doing “more” calculations. The bulk of that is likely the small dimension `k=10` giving a meaningful opportunity to beat BLAS.

EDIT: I think there may also be an indexing bug in the expression above. It looks like something more like

```julia
D[i] += A[j,i]' * B[j,k] * A[k,i]

```

is correct but I derived this using different variables and might have messed up the translation.

---

_[View the full topic](https://discourse.julialang.org/t/fast-diag-a-b-a/98216)._