# Extending broadcasts to matrix exponentials?

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 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> 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> 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
[2] _bcs (repeats 2 times)
[4] combine_axes
[5] instantiate
[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)?

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

1 Like

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

``````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
[2] _bcs (repeats 2 times)
[4] combine_axes
[5] instantiate
[7] top-level scope
@ REPL[107]:1

``````

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> [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> 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.

4 Likes

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

``````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]
``````
1 Like

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

``````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> 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> (cs[:])[3]
2×2 Matrix{Float64}:
2.5   0.44
0.66  1.62
``````

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

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

which naturally doesn’t work. I guess you want

``````getindex.(cs, 3)
``````

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.

1 Like

Thanks, that works great.

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