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

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

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.

3 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]