Extending broadcasts to matrix exponentials?

General Usage
, ,
1

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
   @ .\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)?

2

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

1 Like
3

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
   @ .\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
4

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
5

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
6

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
7

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
8

Thanks, that works great.

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