# Fastest way to perform A \* B \* A’

**URL:** <https://discourse.julialang.org/t/fastest-way-to-perform-a-b-a/117386>\
**Category:** General Usage\
**Tags:** question, blas, mkl, linearalgebra, optimization\
**Created:** [July 23, 2024, 1:10pm UTC](https://discourse.julialang.org/t/fastest-way-to-perform-a-b-a/117386 "2024-07-23T13:10:29Z")\
**Posts on this page:** 1\
**Showing post:** 3

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 23, 2024, 1:32pm UTC](https://discourse.julialang.org/t/fastest-way-to-perform-a-b-a/117386/3 "2024-07-23T13:32:59Z")

</div>

> [@lrnv](#):
>
> I would write the loops myself to compute only that subbloc and try to `@turbo` them.

I wouldn’t recommend this. A fast matrix multiply is about a lot more than SIMD-izing the loops — `@turbo` (from LoopVectorization.jl), or using Tullio.jl or similar, won’t be able to do the higher-level transformations (blocking etc.) that are required to get good performance with such large matrices (where the main concern is memory access).

See also: [LoopVec, Tullio losing to Matrix multiplication](https://discourse.julialang.org/t/loopvec-tullio-losing-to-matrix-multiplication/115798)

> [@eveningsilverfox](#):
>
> Eventually, I need only a certain small sub-block of the result `R`, for instance, `R[10:11,12:13]`.

If you need the subblock `R[a, b]` where `a` and `b` are ranges, just do:

```julia
Rab = @views A[a,:] * B * A[b,:]'

```

and then it will use the fast BLAS matrix-multiplication implementation to compute just the subblock that you want.

if the subblock might not be square, then you should do the “small” dimension first:

```julia
Rab = @views length(a) < length(b) ? (A[a,:] * B) * A[b,:]' : A[a,:] * (B * A[b,:]')

```

> [@eveningsilverfox](#):
>
> `A=inv(I-X)`, where `I` is the identity matrix and `X` is another dense matrix.

You can save additional time by only computing the LU factorization `lu(I - X)` rather than the whole matrix inverse, and then use this to compute only the desired subsets of the columns `Aa = A[a,:]` and `Ab = A[b,:]`:

```julia
LU = lu(I - X)

Aa = zeros(eltype(LU), size(X,1), length(a))
for i = 1:length(a); Aa[a[i], i] = 1; end # Aa = I[:, a]
ldiv!(LU, Aa) # Aa = inv(I - X)[:, a]

Ab = zeros(eltype(LU), size(X,1), length(b))
for i = 1:length(b); Ab[b[i], i] = 1; end # Ab = I[:, b]
ldiv!(LU, Ab) # Ab = inv(I - X)[:, b]

Rab = length(a) < length(b) ? (Aa * B) * Ab' : Aa * (B * Ab') # R[a, b]

```

(Note that you can do even better if `a` and `b` overlap, in which case you can re-use the overlapping portion of the `Aa` columns in `Ab`. But the dominant cost will be the `lu` call if `a` and `b` are small, so re-computing the overlap won’t matter.)

See also: [Why re-use factorization(A) rather than inv(A) - #3 by stevengj](https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234/3) : _95% of the time, if you are computing an explicit matrix inverse, you are doing something suboptimal._

---

_[View the full topic](https://discourse.julialang.org/t/fastest-way-to-perform-a-b-a/117386)._