# Best way to calculate B' \* A \* B

**URL:** <https://discourse.julialang.org/t/best-way-to-calculate-b-a-b/118382>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [August 19, 2024, 6:17pm UTC](https://discourse.julialang.org/t/best-way-to-calculate-b-a-b/118382 "2024-08-19T18:17:04Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [August 19, 2024, 6:17pm UTC](https://discourse.julialang.org/t/best-way-to-calculate-b-a-b/118382/1 "2024-08-19T18:17:04Z")

</div>

I’d like to calculate `B' * A * B`, where `A` is `Hermitian`, without allocating intermediary results. Note that `B` will be, in general, a rectangular matrix. I’m currently using [Tullio](https://github.com/mcabbott/Tullio.jl) to do the calculation, but I was wondering if there are any better ways. Thanks in advance!

---

<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:** [August 19, 2024, 7:09pm UTC](https://discourse.julialang.org/t/best-way-to-calculate-b-a-b/118382/2 "2024-08-19T19:09:22Z")

</div>

Surely you must have at least one allocation for the result, unless you have that pre-allocated. In either case, why not allocate the intermediate result `A * B` as well? Why is it imperative this happen without allocation? If you allocate it outside of any loops, it’s a one-time cost that won’t cause any noticeable slowdown (except at extremely small matrix sizes).

I don’t think there’s an efficient way to do this without additional storage. A non-allocating matrix-matrix-matrix multiply is normally a O(N^4) operation, which is much worse than the O(N^3) allocating version of performing it as a sequence of matrix-matrix multiplies.

You might do slightly better if you have a factorization of `A`, especially a symmetric factorization like Cholesky. If `A == U' * U` then `B' * A * B == (U * B)' * (U * B)`. `lmul!(U, B)` can overwrite `B` in-place with the result of `U * B` for some types of `U`. But this still would likely not save you much (and any gains would be overwhelmed by the cost of factorizing `A` if it was not already). And you’d still be left needing to allocate space for the final result.