# Fast square root of 3 x 3 matrix

**URL:** <https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378>\
**Category:** General Usage\
**Tags:** linearalgebra, staticarrays\
**Created:** [October 27, 2022, 2:57pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378 "2022-10-27T14:57:53Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![AnshPatel7](https://avatars.discourse-cdn.com/v4/letter/a/a698b9/32.png) [@AnshPatel7](https://discourse.julialang.org/u/AnshPatel7)\
**Post date:** [October 27, 2022, 2:57pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/1 "2022-10-27T14:57:53Z")

</div>

Hi,

I’m currently using `sqrt()` from `LinearAlgebra.jl` to compute the square-root of 3 x 3 matrices repeatedly. But this has created a performance bottleneck. So I was wondering if there an implementation of a matrix square root that specializes for 3 x 3 matrices like the ones in `StaticArrays.jl`(which are limited to 2 x 2 matrices). FYI I’m dealing with real, symmetric and positive-definite matrices.

Thanks!

---

<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:** [October 27, 2022, 3:24pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/2 "2022-10-27T15:24:19Z")

</div>

I think the best result will be building your own from the eigen decomposition which StaticArrays has a specialized method for 3x3

---

<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:** [October 27, 2022, 3:29pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/3 "2022-10-27T15:29:41Z")

</div>

Depending on exactly what you need, you might be able to use the Cholesky decomposition, which computes a lower triangular matrix `L` or upper triangular matrix `U` such that `L*L' == U'*U == A`. I believe `StaticArrays` has specialized implementations for `cholesky(A::SMatrix)`. This is similar to a sqrt in some ways, although not in the sense that `S = sqrt(A)` such that `S*S == A` (no transpose needed here and `S` is PSD if `A` is).

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [October 27, 2022, 4:13pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/4 "2022-10-27T16:13:29Z")

</div>

You will get some StaticArrays speedup by default, as long as you remember to wrap your array in a `Symmetric`:

```julia
r_ = rand(3, 3)
r = r_ + r_'

R = r^2
Rh = Symmetric(R)
Rs = SMatrix{3,3}(R)
Rsh = Symmetric(Rs)

@btime sqrt($R); # normal Matrix
  3.538 μs (14 allocations: 2.17 KiB)

@btime sqrt($Rh); # Symmetric Matrix
  3.413 μs (14 allocations: 2.17 KiB)

@btime sqrt($Rs); # SMatrix 
  3.487 μs (16 allocations: 2.31 KiB)

@btime sqrt($Rsh); # Symmetric SMatrix
  726.316 ns (0 allocations: 0 bytes)

```

So you don’t get any benefit of StaticArrays unless you also specify `Symmetric` as well.

---

<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:** [October 27, 2022, 5:00pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/5 "2022-10-27T17:00:25Z")

</div>

> [@mikmoore](#):
>
> Depending on exactly what you need, you might be able to use the Cholesky decomposition

Right, it would be good to know the context here. In many problems where you conceptually might want a matrix square root, there is a way to reformulate it to use Cholesky, which is way cheaper to compute.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [October 27, 2022, 7:08pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/6 "2022-10-27T19:08:04Z")

</div>

> [@Oscar\_Smith](#):
>
> I think the best result will be building your own from the eigen decomposition which StaticArrays has a specialized method for 3x3

The `Symmetric` LinearAlgebra `sqrt` is implemented with `eigen`, so with a `Symmetric SMatrix` you don’t have to build it yourself.

If speed is more critical than precision you can do considerably better than that, however. This function only works for sufficiently distinct eigenvalues, so needs to be complemented with some special cases for colliding eigenvalues.

```julia
function f(x)
    l = eigvals(x)
    J = Symmetric(SMatrix{3, 3}(I))
    A1 = Symmetric(x - l[1] * J)
    A2 = Symmetric(x - l[2] * J)
    A3 = Symmetric(x - l[3] * J)
    B3 = Symmetric(A1 * A2) / ((l[3] - l[1]) * (l[3] - l[2]))
    B2 = Symmetric(A1 * A3) / ((l[2] - l[1]) * (l[2] - l[3]))
    B1 = Symmetric(A2 * A3) / ((l[1] - l[2]) * (l[1] - l[3]))
    return sqrt(l[1]) * B1 + sqrt(l[2]) * B2 + sqrt(l[3]) * B3
end

```

With a `Symmetric{Smatrix}` I measure it 5-6 times faster than `sqrt` and it probably leaves a fair amount of micro-optimizations on the table.

---

<div class="post-metadata">

**Author:** ![AnshPatel7](https://avatars.discourse-cdn.com/v4/letter/a/a698b9/32.png) [@AnshPatel7](https://discourse.julialang.org/u/AnshPatel7)\
**Post date:** [October 28, 2022, 1:47pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/7 "2022-10-28T13:47:55Z")

</div>

Thanks this works, even when R is a `MMatrix{3, 3}`.

---

<div class="post-metadata">

**Author:** ![AnshPatel7](https://avatars.discourse-cdn.com/v4/letter/a/a698b9/32.png) [@AnshPatel7](https://discourse.julialang.org/u/AnshPatel7)\
**Post date:** [October 28, 2022, 1:58pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/8 "2022-10-28T13:58:21Z")

</div>

I need to factorize a 3 x 3 matrix A such that B\*B’ = A and then B is used elsewhere in the code. The problem with the Cholesky decomposition is that now I’m not sure that my A will always be positive definite. But it is true that Cholesky is much faster than computing the square-root even with StaticArrays and Symmetric.

---

<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:** [October 28, 2022, 2:05pm UTC](https://discourse.julialang.org/t/fast-square-root-of-3-x-3-matrix/89378/9 "2022-10-28T14:05:04Z")

</div>

Cholesky is definitely what you want here.