# Most efficient way of computing B'\* inv(A)\* B for A large

**URL:** <https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513>\
**Category:** Specific Domains\
**Tags:** linearalgebra\
**Created:** [September 29, 2020, 9:31pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513 "2020-09-29T21:31:28Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 9:31pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/1 "2020-09-29T21:31:28Z")

</div>

I have a positive definite high-dimensional matrix A and wish to compute B’ \* inv(A) \* B where the result is low-dimensional. What is the fastest way of accomplishing this for generic B?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 29, 2020, 9:33pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/2 "2020-09-29T21:33:33Z")

</div>

To start, writing it as `(A\B')*B` will be more accurate and likely faster.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 29, 2020, 9:52pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/3 "2020-09-29T21:52:11Z")

</div>

Do this: `B' * (A \ B)`

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 9:59pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/4 "2020-09-29T21:59:11Z")

</div>

Sure. I was looking for something a little better than that. And preferably something that guarantees that the result is exactly symmetric (instead of symmetric up to rounding error).

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 29, 2020, 10:04pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/5 "2020-09-29T22:04:45Z")

</div>

`Symmetric(B' * (A \ B))` would do that.

---

<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:** [September 29, 2020, 10:12pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/6 "2020-09-29T22:12:41Z")

</div>

Is `A` hermitian?

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 10:16pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/7 "2020-09-29T22:16:35Z")

</div>

It’s positive-definite even. Are you thinking a Cholesky decomposition?

---

<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:** [September 29, 2020, 10:16pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/8 "2020-09-29T22:16:50Z")

</div>

Yes.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 10:18pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/9 "2020-09-29T22:18:47Z")

</div>

That would work. I’d been hoping for something canned in MKL or some such library for the whole thing.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 10:21pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/10 "2020-09-29T22:21:48Z")

</div>

I’ll see if I can do Cholesky, then Tullio.

---

<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:** [September 29, 2020, 10:26pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/11 "2020-09-29T22:26:07Z")

</div>

Yeah, Cholesky and then Tullio was exactly my thought. 🙂

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 29, 2020, 10:48pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/12 "2020-09-29T22:48:24Z")

</div>

Is `B` sparse? Is `A` sparse?

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 10:56pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/13 "2020-09-29T22:56:26Z")

</div>

No no

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [September 29, 2020, 11:08pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/14 "2020-09-29T23:08:20Z")

</div>

A = LL^T, C = L^{-1} B, then compute C^T C.

Edit: by L^{-1} B I mean `L \ B`.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 11:09pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/15 "2020-09-29T23:09:36Z")

</div>

Right, that’s Cholesky.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [September 29, 2020, 11:11pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/16 "2020-09-29T23:11:27Z")

</div>

It’s Cholesky, but don’t do `A \ B`, it’s unnecessary.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 29, 2020, 11:12pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/17 "2020-09-29T23:12:07Z")

</div>

Isn’t that `L \ B`?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [September 29, 2020, 11:35pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/18 "2020-09-29T23:35:01Z")

</div>

I see: you are saying that only backward substitution is needed as opposed to back/forward for `A \ B`.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [September 29, 2020, 11:46pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/19 "2020-09-29T23:46:08Z")

</div>

Yes, exactly.

---

<div class="post-metadata">

**Author:** ![Joris\_Pinkse](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joris_pinkse/32/216398_2.png) [@Joris\_Pinkse](https://discourse.julialang.org/u/Joris_Pinkse)\
**Post date:** [September 29, 2020, 11:56pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513/20 "2020-09-29T23:56:54Z")

</div>

For a 2,000 by 2,000 matrix A and 2,000 by 10 matrix B, I tried

```julia
function one()
    C[:,:] =Symmetric( B' * (A \ B ) )
end

function two()
    cholesky!(A)
    D[:,:] = A \ B
    @tullio C[i,j] = D[k,i] * D[k,j]
end

function three()
    D[:,:] = cholesky(A).U \ B
    @tullio C[i,j] = D[k,i] * D[k,j]
end

function four()
    D[:,:] = cholesky(A).U \ B
    C[:,:] = D' * D
end
 
function five()
    D = cholesky(A).U \ B
    C[:,:] = D' * D
end

```

Here are the second run times:  
0.032261 seconds (12 allocations: 30.687 MiB, 8.10% gc time)  
0.031080 seconds (25 allocations: 30.686 MiB)  
0.015401 seconds (24 allocations: 30.671 MiB, 13.39% gc time)  
0.009717 seconds (15 allocations: 30.671 MiB)  
0.008867 seconds (13 allocations: 30.671 MiB)

I’m puzzled by the garbage collection time for three(), but there appears to be a clear winner. My matrices are of the order 4\_000 \* 4\_000 and it won’t have to be done more than a few thousand times per run, so this will be good enough.

Thanks all!  
(there were some preallocations there)

[Next page](https://discourse.julialang.org/t/most-efficient-way-of-computing-b-inv-a-b-for-a-large/47513.md?page=2)