How to add other arguments to function turned to LinearOperator?

[erny123]

Hi all,

Trying to write a matrix-free solver using LinearOperators.jl and Krylov.jl
However, from this: https://jso.dev/tutorials/introduction-to-linear-operators/
I'm seeing that you can only use functions that follow the 5-argument (and 3) mul! operation.

But I would like to create a linear operator that is able to take in more variables such as:

function customfunc!(out, in, a, b, dt, nu)
for i in 1:length(in)
out[i] = in[i]*dt*(nu[i]+nu[i+1])/2.0
end

In LinearMaps.jl it would be something like:

customfunc = (N, dt, nu) -> LinearMap(N; ismutating=true) do Z,U

for i in 1:N
Z[i] = U[i]*dt*(nu[i]+nu[i+1])/2.0
end
return Z
end

CFunc = customfunc(100, 0.01, [1.0,2.0,3.0])

Then I could use CFunc in IterativeSolvers.jl.

Any suggestions?

[geoffroyleconte]

Hi, do you mean something like this:

  K = LinearOperator(
    T,
    id.nvar + id.ncon,
    id.nvar + id.ncon,
    true,
    true,
    (res, v, α, β) ->
      opK2_5prod!(res, id.nvar, fd.Q, D, fd.A, sqrtX1X2, tmp1, tmp2, δv, v, α, β, fd.uplo),
  )

(taken from https://github.com/JuliaSmoothOptimizers/RipQP.jl/blob/a44357fa457a0d6a38337e18c49fd4b72fdab3b5/src/iterations/solvers/augmented/K2_5Krylov.jl#L159)?
Here if think opK2_5prod! is your customfunc!.

If you want to modify the values of the vector sqrtX1X2 for example, it will be modified in the function opK2_5prod!.
If you want to do this with scalars it is not possible, but you can use for example δv = [1.0] and modify δv[1] if you need.
If δv stays always the same for the operator K, you can use an argument with a default value.

Another option (maybe more advanced) is to write a subtype of AbstractLinearOperator{T}.
See for example:

LinearOperators.jl/src/DiagonalHessianApproximation.jl

Lines 17 to 52 in 9cb279b

	mutable struct DiagonalQN{T <: Real, I <: Integer, V <: AbstractVector{T}, F} <:
	AbstractDiagonalQuasiNewtonOperator{T}
	d::V # Diagonal of the operator
	nrow::I
	ncol::I
	symmetric::Bool
	hermitian::Bool
	prod!::F
	tprod!::F
	ctprod!::F
	nprod::I
	ntprod::I
	nctprod::I
	args5::Bool
	use_prod5!::Bool # true for 5-args mul! and for composite operators created with operators that use the 3-args mul!
	allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
	psb::Bool
	end
	"""
	DiagonalQN(d)
	Construct a linear operator that represents a diagonal quasi-Newton approximation.
	The approximation satisfies the weak secant equation and is not guaranteed to be
	positive definite.
	# Arguments
	- `d::AbstractVector`: initial diagonal approximation;
	- `psb::Bool`: whether to use the diagonal PSB update or the Andrei update.
	"""
	function DiagonalQN(d::AbstractVector{T}, psb::Bool = false) where {T <: Real}
	prod = (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β)
	n = length(d)
	DiagonalQN(d, n, n, true, true, prod, prod, prod, 0, 0, 0, true, true, true, psb)
	end

[erny123]

@geoffroyleconte so from your reply it's not exactly what I need although the custom AbstractLinearOperator might be worth looking into.

I ended up doing the following which I posted on the Julia discord: https://discourse.julialang.org/t/how-to-add-other-arguments-to-function-turned-to-linearoperator/110475/4

You can use an anonymous function inside another anonymous function such as:

customfunc = (N, dt, nu) -> (Z,U,alpha,beta) -> begin

for i in 1:N
Z[i] = U[i]*dt*(nu[i]+nu[i+1])/2.0
end
return Z
end
Then you initialize with your parameters:

N =12
dt = 0.5
nu=0.01
custmul! = customfun(N,dt,nu)

Then you can put this into a linear operator because this is now a 4 argument function that works as:

custmul!(Out, In, alpha, beta)

Do you think this is a somewhat optimal way to do it? It might also be a quick lazy way

[geoffroyleconte]

Great if this works!

I was thinking about something like

function customfunc!(out, in, a, b, dt, nu)
  for i in 1:length(in)
  out[i] = in[i]*dt*(nu[i]+nu[i+1])/2.0
  end
  out
end
dt = 0.5
nu = [1.0,2.0,3.0]
custmul! = (out, in, alpha,beta) -> customfunc!(out, in, a, b, dt, nu)

or

dt = 0.5
nu = [1.0,2.0,3.0]
function custmul!(out, in, a, b; dt = dt, nu = nu)
  for i in 1:length(in)
  out[i] = in[i]*dt*(nu[i]+nu[i+1])/2.0
  end
  out
end

but I think this is similar to your answer.