# Matrix power not memory optimal

**URL:** <https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201>\
**Category:** Profiling\
**Tags:** performance, linearalgebra, memory-allocation\
**Created:** [June 25, 2024, 12:46pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201 "2024-06-25T12:46:48Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [June 25, 2024, 12:46pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/1 "2024-06-25T12:46:48Z")

</div>

It seems that computing A^p with integer p for a dense matrix A with `^(A,p)` is not memory allocation optimal:

```julia-repl
julia> using BenchmarkTools, LinearAlgebra
julia> A=randn(256,256); A=A/norm(A);
julia> @btime A8=A^8;
  813.054 μs (6 allocations: 1.50 MiB)
julia> @btime begin; 
    A2=similar(A); mul!(A2,A,A); 
    A4=similar(A); mul!(A4,A2,A2); 
    A8=A2; mul!(A8,A4,A4); end;
  729.419 μs (4 allocations: 1.00 MiB)
julia> @btime ((A^2)^2)^2;
  752.717 μs (6 allocations: 1.50 MiB)

```

In the second btime, I recycle the memory slots when repeating the powers. The difference should be even larger for higher powers. As far as I can tell, the reason is that `^(A,::Integer)` makes a call to `Base.power_by_squaring` in `intfuncs.jl` which does not take such memory aspects into account (in particular line 330)

> <https://github.com/JuliaLang/julia/blob/5654e6043823717e085239f6509413410106e902/base/intfuncs.jl#L314-L334>

Is there a good repeated squaring in some package in Julia?

---

<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:** [June 25, 2024, 3:25pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/2 "2024-06-25T15:25:16Z")

</div>

this would be a pretty good first pr for someone. the fix is quite straightforward

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [June 25, 2024, 3:53pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/3 "2024-06-25T15:53:39Z")

</div>

The function `power_by_squaring` seems to be written for scalars and immutable types, I don’t see an easy fix where `mul` is replaced by `mul!`. One can make some conditions on the type, but since this is in `Base` references to any `Matrix` is not ideal. What fix did you have in mind?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 25, 2024, 3:55pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/4 "2024-06-25T15:55:06Z")

</div>

> [@Oscar\_Smith](#):
>
> this would be a pretty good first pr for someone. the fix is quite straightforward

How would it handle immutable matrices like FillArrays.jl or StaticArrays.jl?

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [June 25, 2024, 3:58pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/5 "2024-06-25T15:58:08Z")

</div>

> [@gdalle](#):
>
> How would it handle immutable matrices like [FillArrays.jl](https://juliahub.com/ui/Packages/General/FillArrays) or [StaticArrays.jl](https://juliahub.com/ui/Packages/General/StaticArrays)?

`StaticArrays` are not necessarily immutable. Not sure we have many immutable matrices, but I get your point.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 25, 2024, 4:01pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/6 "2024-06-25T16:01:53Z")

</div>

> [@jarl](#):
>
> Not sure we have many immutable matrices, but I get your point.

They hide in many places

```julia
julia> using LinearAlgebra

julia> I
UniformScaling{Bool}
true*I

julia> (2I)^3
UniformScaling{Int64}
8*I

```

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [June 25, 2024, 4:09pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/7 "2024-06-25T16:09:18Z")

</div>

Since `I` is not a `AbstractMatrix`, your example does not speak against a type-check. I believe your call reduces to a call to `power_by_squaring(::Float64,::Int)`, which should remain fine.

Not saying I have a good solution though…

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 25, 2024, 4:11pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/8 "2024-06-25T16:11:34Z")

</div>

> [@jarl](#):
>
> Since `I` is not a `AbstractMatrix`, your example does not speak against a type-check.

My bad. But I’m not sure it’s easy to statically deduce whether a matrix is mutable or not?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 25, 2024, 4:48pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/9 "2024-06-25T16:48:31Z")

</div>

```julia
julia> using StaticArrays
WARNING: using StaticArrays.pop in module Main conflicts with an existing identifier.

julia> A = @SMatrix rand(3,3);

julia> B = similar(A); B .= 3;

julia> convert(typeof(parent(A)), B)
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 3.0 3.0 3.0
 3.0 3.0 3.0
 3.0 3.0 3.0

```

Can use `similar`, mutate it, and then use `convert` to return an answer of the correct type.

`typeof(parent(A))` is so that this would work for at least some wrappers like `view`.  
Alternatively, could do a more general

```julia
RT = Base.promote_op() do
  # capture `A`, not a type!
  convert(typeof(parent(A)), B)
end
if RT === typeof(parent(A))
    # only try converting if it is known to return the expected type at compile time
    convert(typeof(parent(A)), B)
else
    B
end

```

An easier approach is probably to introduce some interface.

---

<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:** [June 25, 2024, 5:27pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/10 "2024-06-25T17:27:47Z")

</div>

> [@Elrod](#):
>
> Can use `similar`, mutate it, and then use `convert` to return an answer of the correct type.

In general, the return type for `^(A::AbstractMatrix,::Integer)` should probably be what is returned by `similar(A)` (which is always mutable)? StaticArrays can overload `^(::SMatrix, ::Integer)` to call `power_by_squaring` directly and avoid allocating a mutable copy.

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [June 25, 2024, 5:45pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/11 "2024-06-25T17:45:00Z")

</div>

The example shows how the performance can be improved for `Matrix`. The same does not hold for `SparseMatrix` because of increased fill-in, not captured well with `similar`.

```julia
julia> A=sprandn(256,256,0.1); A=A/norm(A);
julia> @btime A8=A^8;
  34.995 ms (19 allocations: 3.01 MiB)
julia> @btime begin; 
          A2=similar(A); mul!(A2,A,A); 
          A4=similar(A); mul!(A4,A2,A2); 
          A8=A2; mul!(A8,A4,A4); end;
  205.170 ms (51 allocations: 3.26 MiB)

```

Therefore, I’m currently I would say it’s better to have a separate implementation specifically for dense matrices, i.e., `^(::Matrix,::Int)`, and keep the current as fallback eg for `SMatrix` and `SparseMatrix`.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [June 25, 2024, 6:12pm UTC](https://discourse.julialang.org/t/matrix-power-not-memory-optimal/116201/12 "2024-06-25T18:12:32Z")

</div>

> [@jarl](#):
>
> ^(::Matrix,::Int)

Perhaps this should be specialized for a `StridedMatrix`, which should cover most dense types.