# Extending broadcasts to matrix exponentials?

**URL:** <https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384>\
**Category:** General Usage\
**Tags:** proposal, broadcast, matrices\
**Created:** [May 2, 2022, 8:04pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384 "2022-05-02T20:04:35Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![g2g](https://avatars.discourse-cdn.com/v4/letter/g/41988e/32.png) [@g2g](https://discourse.julialang.org/u/g2g)\
**Post date:** [May 2, 2022, 8:04pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/1 "2022-05-02T20:04:35Z")

</div>

I’m fairly new to Julia and I have not yet looked at Julia’s source code. I’m wondering

- whether extending broadcasts to approximate [matrix exponentials](https://en.wikipedia.org/wiki/Matrix_exponential) would be an improvement to the language, and
- whether the code change would not be too difficult to implement.

Currently (in version 1.7.2),

```julia
julia> B = [1 1/factorial(1) 1/factorial(2) 1/factorial(3) 1/factorial(4) 1/factorial(5) 1/factorial(6)]
1×7 Matrix{Float64}:
 1.0 1.0 0.5 0.166667 0.0416667 0.00833333 0.00138889

julia> x = 0.2
0.2

julia> A = [x^0; x; x^2; x^3; x^4; x^5; x^6]
7-element Vector{Float64}:
 1.0
 0.2
 0.04000000000000001
 0.008000000000000002
 0.0016000000000000003
 0.0003200000000000001
 6.400000000000002e-5

julia> c = B * A
1-element Vector{Float64}:
 1.2214027555555556

julia> exp(0.2)
1.2214027581601699

julia> c = B .* A
7×7 Matrix{Float64}:
 1.0 1.0 0.5 0.166667 0.0416667 0.00833333 0.00138889
 0.2 0.2 0.1 0.0333333 0.00833333 0.00166667 0.000277778
 0.04 0.04 0.02 0.00666667 0.00166667 0.000333333 5.55556e-5
 0.008 0.008 0.004 0.00133333 0.000333333 6.66667e-5 1.11111e-5
 0.0016 0.0016 0.0008 0.000266667 6.66667e-5 1.33333e-5 2.22222e-6
 0.00032 0.00032 0.00016 5.33333e-5 1.33333e-5 2.66667e-6 4.44444e-7
 6.4e-5 6.4e-5 3.2e-5 1.06667e-5 2.66667e-6 5.33333e-7 8.88889e-8

```

So when x = 0.2, B \* A gives an approximation of exp(0.2). But when x is a square matrix such as [1 4; 1 1], B \* A and B .\* A result in errors, not an approximation of exp([1 4; 1 1]).

```julia
julia> x = [1 4; 1 1]
2×2 Matrix{Int64}:
 1 4
 1 1

julia> exp(x)
2×2 Matrix{Float64}:
 10.2267 19.7177
  4.92941 10.2267

julia> A = [x^0; x; x^2; x^3; x^4; x^5; x^6]
14×2 Matrix{Int64}:
   1 0
   0 1
   1 4
   1 1
   5 8
   2 5
  13 28
   7 13
  41 80
  20 41
 121 244
  61 121
 365 728
 182 365

julia> C = B .* A
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 7 and 2")
Stacktrace:
 [1] _bcs1
   @ .\broadcast.jl:516 [inlined]
 [2] _bcs (repeats 2 times)
   @ .\broadcast.jl:510 [inlined]
 [3] broadcast_shape
   @ .\broadcast.jl:504 [inlined]
 [4] combine_axes
   @ .\broadcast.jl:499 [inlined]
 [5] instantiate
   @ .\broadcast.jl:281 [inlined]
 [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(*), Tuple{Matrix{Float64}, Matrix{Int64}}})
   @ Base.Broadcast .\broadcast.jl:860
 [7] top-level scope
   @ REPL[105]:1

```

Wouldn’t it be more consistent if the elements of A were recognized as square matrices, each to be broadcast multiplied by the corresponding element in B, so that B .\* A is again an approximation of exp(x)?

---

<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:** [May 2, 2022, 8:08pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/2 "2022-05-02T20:08:17Z")

</div>

this works, you just need to not use ; to construct the array

---

<div class="post-metadata">

**Author:** ![g2g](https://avatars.discourse-cdn.com/v4/letter/g/41988e/32.png) [@g2g](https://discourse.julialang.org/u/g2g)\
**Post date:** [May 2, 2022, 8:16pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/3 "2022-05-02T20:16:03Z")

</div>

I don’t quite understand your suggestion. If I instead try to construct A using a transpose, I get the same error.

```julia
julia> A = [x^0 x x^2 x^3 x^4 x^5 x^6]'
14×2 adjoint(::Matrix{Int64}) with eltype Int64:
   1 0
   0 1
   1 1
   4 1
   5 2
   8 5
  13 7
  28 13
  41 20
  80 41
 121 61
 244 121
 365 182
 728 365

julia> C = B .* A
ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 7 and 2")
Stacktrace:
 [1] _bcs1
   @ .\broadcast.jl:516 [inlined]
 [2] _bcs (repeats 2 times)
   @ .\broadcast.jl:510 [inlined]
 [3] broadcast_shape
   @ .\broadcast.jl:504 [inlined]
 [4] combine_axes
   @ .\broadcast.jl:499 [inlined]
 [5] instantiate
   @ .\broadcast.jl:281 [inlined]
 [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(*), Tuple{Matrix{Float64}, LinearAlgebra.Adjoint{Int64, Matrix{Int64}}}})
   @ Base.Broadcast .\broadcast.jl:860
 [7] top-level scope
   @ REPL[107]:1

```

---

<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:** [May 2, 2022, 8:23pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/4 "2022-05-02T20:23:47Z")

</div>

Both `[a; b; c]` and `[a b c]` do _concatenation_, the former vertical, the latter horizontal. You don’t want concatenation here, you want to build a vector of matrices. The correct way to build a vector is `[a, b, c]` (note: comma, not semicolon or space):

```julia
julia> [x^0, x, x^2, x^3, x^4, x^5, x^6]
7-element Vector{Matrix{Int64}}:
 [1 0; 0 1]
 [1 4; 1 1]
 [5 8; 2 5]
 [13 28; 7 13]
 [41 80; 20 41]
 [121 244; 61 121]
 [365 728; 182 365]

```

You can also use broadcasting:

```julia
julia> Ref(x).^(0:6)
7-element Vector{Matrix{Int64}}:
 [1 0; 0 1]
 [1 4; 1 1]
 [5 8; 2 5]
 [13 28; 7 13]
 [41 80; 20 41]
 [121 244; 61 121]
 [365 728; 182 365]

```

The `Ref` around `x` is necessary to ensure that it is treated like a scalar, and not broadcast over the individual elements.

---

<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:** [May 2, 2022, 8:27pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/5 "2022-05-02T20:27:06Z")

</div>

Here’s the difference between concatenation and vector literal construction:

```julia
julia> x = [1,2,3];

julia> y = [4,5,6];

julia> [x; y]
6-element Vector{Int64}:
 1
 2
 3
 4
 5
 6

julia> [x y]
3×2 Matrix{Int64}:
 1 4
 2 5
 3 6

julia> [x, y]
2-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [4, 5, 6]

```

---

<div class="post-metadata">

**Author:** ![g2g](https://avatars.discourse-cdn.com/v4/letter/g/41988e/32.png) [@g2g](https://discourse.julialang.org/u/g2g)\
**Post date:** [June 8, 2022, 9:02pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/6 "2022-06-08T21:02:16Z")

</div>

> [@DNF](#):
>
> You can also use broadcasting:
> 
> `julia> Ref(x).^(0:6)`

Using your broadcasting tip, I was glad to see that cumsum() can be used to calculate the matrix geometric series partial sums:

```julia
julia> x = [0.8 0.2; 0.3 0.4]
2×2 Matrix{Float64}:
 0.8 0.2
 0.3 0.4

julia> s = Ref(x).^(0:10)
11-element Vector{Matrix{Float64}}:
 [1.0 0.0; 0.0 1.0]
 [0.8 0.2; 0.3 0.4]
 [0.7000000000000002 0.24000000000000005; 0.36 0.22000000000000003]
 [0.6320000000000001 0.23600000000000007; 0.35400000000000004 0.16000000000000003]
 [0.5764000000000002 0.22080000000000008; 0.33120000000000005 0.13480000000000003]
 [0.5273600000000002 0.2036000000000001; 0.3054000000000001 0.12016000000000004]
 [0.4829680000000003 0.1869120000000001; 0.2803680000000001 0.10914400000000005]
 [0.4424480000000003 0.17135840000000008; 0.25703760000000014 0.09973120000000005]
 [0.4053659200000003 0.1570329600000001; 0.2355494400000001 0.09130000000000005]
 [0.37140262400000035 0.1438863680000001; 0.21582955200000015 0.08362988800000004]
 [0.3402880096000004 0.1318350720000001; 0.19775260800000014 0.07661786560000006]

julia> cs = cumsum(s, dims=1)
11-element Vector{Matrix{Float64}}:
 [1.0 0.0; 0.0 1.0]
 [1.8 0.2; 0.3 1.4]
 [2.5 0.44000000000000006; 0.6599999999999999 1.62]
 [3.1320000000000006 0.6760000000000002; 1.014 1.7800000000000002]
 [3.708400000000001 0.8968000000000003; 1.3452000000000002 1.9148]
 [4.235760000000001 1.1004000000000003; 1.6506000000000003 2.0349600000000003]
 [4.718728000000001 1.2873120000000005; 1.9309680000000005 2.1441040000000005]
 [5.161176000000002 1.4586704000000006; 2.1880056000000008 2.2438352000000004]
 [5.566541920000002 1.6157033600000008; 2.423555040000001 2.3351352000000007]
 [5.937944544000002 1.7595897280000008; 2.639384592000001 2.4187650880000007]
 [6.278232553600002 1.891424800000001; 2.8371372000000012 2.495382953600001]

```

Is there some indexing tip to slice out a particular element of each matrix (e.g., to see how the element is converging)? Extracting one matrix and getting an element works fine:

```julia
julia> cs[2]
2×2 Matrix{Float64}:
 1.8 0.2
 0.3 1.4

julia> (cs[2])[3]
0.2

```

However, my attempt at slicing an element from each matrix is incorrect:

```julia
julia> (cs[:])[3]
2×2 Matrix{Float64}:
 2.5 0.44
 0.66 1.62

```

---

<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:** [June 8, 2022, 10:41pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/7 "2022-06-08T22:41:01Z")

</div>

> [@g2g](#):
>
> `julia> (cs[:])[3]`

`cs[:]` just creates a copy of `cs`, so the above is no different from writing

```julia
cs2 = copy(cs)
cs2[3]

```

which naturally doesn’t work. I guess you want

```julia
getindex.(cs, 3)

```

> [@g2g](#):
>
> `cs = cumsum(s, dims=1)`

Since `s` is a vector, you don’t need the `dims` keyword here, just write `cs = cumsum(s)`.

But if you don’t want to waste too much performance, I wouldn’t recommend broadcasting the exponential, since that is extremely wasteful (each iteration takes longer than the previous). Instead use a loop, and produce the next element using a single multiplication with the previous.

---

<div class="post-metadata">

**Author:** ![g2g](https://avatars.discourse-cdn.com/v4/letter/g/41988e/32.png) [@g2g](https://discourse.julialang.org/u/g2g)\
**Post date:** [June 8, 2022, 11:18pm UTC](https://discourse.julialang.org/t/extending-broadcasts-to-matrix-exponentials/80384/8 "2022-06-08T23:18:07Z")

</div>

Thanks, that works great.

From my beginner perspective, one of my initial attempts was cs.[3] but that caused an invalid syntax error.