Parametrized maps

[andreasvarga]

Is any way to define linear maps depending on some parameters? For example, I would like to define a time-dependent map A(t) or even a map A(t,p) depending on time values t and some parameters inp, to be used in DifferentialEquations.jl to define the function as f(u,p,t) = A(t,p)*u ?

[dkarrasch]

DifferentialEquations.jl has lots of tooling around defining (possibly time- and state-dependent) linear operators, though I'm not sure if they always need to be backed by a matrix. In which way is your operator defined? As a function that takes parameters t and p? Then you could write a function A(t, p) = LinearMap(myfun(t, p), m, n; properties...). Could you share a few more details on the context?

[andreasvarga]

Sorry if I am wrong, but here is what I am expecting. For a matrix based linear map I have:

julia> B = LinearMap(rand(2,2))
2×2 LinearMaps.WrappedMap{Float64} of
  2×2 Matrix{Float64}

julia> B*[1;1]
2-element Vector{Float64}:
 0.461044462571232
 1.3866005501303436

For a time-dependent linear map, I tried:

julia> a(t) = [cos(t) sin(t); 1 -1]
a (generic function with 1 method)

julia> A(t) = LinearMap(a(t), 2, 2; ismutating = false)
A (generic function with 1 method)

julia> A(1)
2×2 LinearMaps.FunctionMap{Float64}([0.5403023058681398 0.8414709848078965; 1.0 -1.0]; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)

julia> A(1)*[1;1]
ERROR: MethodError: objects of type Matrix{Float64} are not callable
Use square brackets [] for indexing an Array.
Stacktrace:
 [1] *(A::LinearMaps.FunctionMap{Float64, Matrix{Float64}, Nothing}, x::Vector{Int64})
   @ LinearMaps C:\Users\Andreas\.julia\packages\LinearMaps\1cWDb\src\functionmap.jl:45
 [2] top-level scope
   @ REPL[42]:1

What I am doing wrong?

Actually

julia> A(1)*B
2×2 LinearMaps.CompositeMap{Float64} with 2 maps:
  2×2 LinearMaps.WrappedMap{Float64} of
    2×2 Matrix{Float64}
  2×2 LinearMaps.FunctionMap{Float64}([0.5403023058681398 0.8414709848078965; 1.0 -1.0]; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)

julia> (A(1)*B)*[1;1]
ERROR: MethodError: objects of type Matrix{Float64} are not callable
Use square brackets [] for indexing an Array.

Regarding your question: I would like to perform composition, linear combination, concatenations, etc. operations on time-dependent operators. This would be helpful in manipulating linear time-varying dynamical systems of the form

dx(t)/dt = A(t)*x(t)+B(t)*u(t)
y(t) = C(t)*x(t) + D(t)*u(t)

I just started a new project. In this moment I not yet decided how to define a system object which contains the four matrices and can be easily used for further manipulations (e.g., series, parallel, diagonal, couplings). Sorry bothering you with my problems.

[dkarrasch]

Sorry bothering you with my problems.

Don't worry, that's what we're all here for.

Since a is a function, the LinearMap constructor thinks this should be a FunctionMap, which is not fully accurate since you're constructing a matrix for given t. Anyway, the FunctionMap requires the function argument to be the one that is performing matvec multiplication, so this works:

using LinearMaps, LinearAlgebra
B = LinearMap(rand(2,2))
a(t) = [cos(t) sin(t); 1 -1]
A(t) = LinearMap((y,x) -> mul!(y, a(t), x), 2, 2; ismutating = true)
A(1)
A(1)*[1;1]
A(1)*B
(A(1)*B)*[1;1]

Regarding your question: I would like to perform composition, linear combination, concatenations, etc. operations on time-dependent operators.

Awesome, this is music to my ears.

I'd still recommend to take a look at https://diffeq.sciml.ai/stable/features/diffeq_operator/ or http://diffeqoperators.sciml.ai/stable/. So, in case you write your own type, you may want to consider subtyping to AbstractDiffEqOperator to hook into the ODE solver suite directly. Currently, it only supports linear combinations and composition, but no concatenation, though.

[andreasvarga]

Actually I already used https://diffeq.sciml.ai/stable/features/diffeq_operator/ to implement the function for a particular integration method (Magnus' mthod). I think it will be possible to implement function matrices, as for example G(t) = [A1(t) B1(t)*C2(t); 0 A2(t)](needed in series coupling of systems) directly from the matrix functions A1(t), B1(t), C2(t) and A2(t), without resorting to linear maps. The question is, which approach is more efficient when performing heavy computations involving these matrices.

I am not sure, but I have the impression that each time-evaluation of the operator A(t), which you defined, creates an instance of the operator. Could this lead to some performance loss (e.g., when integrating ODE's)?

[dkarrasch]

I am not sure, but I have the impression that each time-evaluation of the operator A(t)creates an instance of the operator. Could this lead to some performance loss (e.g., when integrating ODE's)?

That's correct. Performancewise this will depend on how much effort it takes to evaluate A(t)*x. If that's rather intense, then creating the thin BlockMap should be neglible. But you're right, at every time evaluation, it will go through the code for size consistency checking etc. I assume that for that reason, DiffEqOperators has this update_coefficient! function to update in-place. Your blocks, are they arrays, or do you work with function-based LinearMaps? In that case, you may want to consider also BlockArrays.jl.

Conversely, if you have rather fixed ways ("patterns") of combining some maps, you could create your own type for those specific combinations. (I'm making something up)

struct SeriesCoupling{T,...} <: LinearMap{T}
    A1::A
    A2::AA
    B1::B
    C2::C
end

and then write specific multiplication routines for it. Or

struct SeriesCoupling{T,As <: NTuple{4,LinearMap} <: LinearMap{T}
    maps::As
end

Depends on how many patterns you will need. Note that matrix-based LinearMaps are just thin wrappers around that matrix. So you can still modify the wrapped matrix elementwise.

[dkarrasch]

There is a new project in development, that does what I was considering to do here at some point: make LinearMaps potentially parameter-dependent, i.e., merge LinearMaps.jl with the usual DiffEq signature (u,p,t), their trait system etc. They also preallocate any caches they may need during the operator application. The package is SciMLOperators.jl. They are working hard to beat LinearMaps.jl performance-wise. 🤣