SciMLOperators “*” function for TransposedOperator and AdjointOperator

[fdekerme]

Hello,
https://discourse.julialang.org/t/scimloperators-function-for-transposedoperator/129422

The dispatch of multiplication (*) for TransposedOperator is missing:

function Afunc!(w,v,u,p,t)
    w[1] = -2v[1] + v[2]
    for i in 2:4
        w[i] = v[i-1] - 2v[i] + v[i+1]
    end
    w[5] = v[4] - 2v[5]
    nothing
end

function Afunc!(v,u,p,t)
    w = zeros(5)
    Afunc!(w,v,u,p,t)
    w
end

mfopA = FunctionOperator(Afunc!, zeros(5), zeros(5))
mfopA_t = transpose(mfopA)
mfopA * rand(5)
mfopA_t * rand(5)`

SciMLOperators v1.3.0
Julia 1.11.5

Thanks

[fdekerme]

Same problem for AdjointOperator

mfopA' * rand(5)

[fdekerme]

Hi. Any news about this issue? I need to transpose an operator and I am stuck with this bug.
Thanks

[ChrisRackauckas]

The transpose operator has the opposite interaction:

rand(5)' * mfopA_t # works!

In order to define the adjoint operator here, you need to do reverse-mode AD of the original function operator, i.e. this is the vjp operation. It's not a trivial thing, though we could use SciMLJacobianOperators.jl (inside NonlinearSolve.jl) in order to generate it.

[fdekerme]

I see, thanks for the answer.
No special case for square matrix/operator ?
For instance, how can we define
$\lambda (W - I)^T(W - I)$ with $W$ and $I$ square matrix?

[ChrisRackauckas]

Why are you doing A'A? Is that for a least squares calculation? If so, you can do the same via a QR factorization or a GMRES directly. There is basically no reason to construct A'A, and it's generally a bad idea because it will square the condition number of the operator.

There isn't a special case when you have a function operator. What we could do is give you a way to supply the transpose operation and use that?

[fdekerme]

I am trying to solve the linear problem $(A^T V^{−1} A+ λ·∇G)x = A^T V^{−1}\mu$ with different regularization $∇G$. I already implemented L2 penalty and Huber loss with SciMLOperators. Now, I want to try the "similarity based" regularization of this paper: https://doi.org/10.1118/1.4947485. I have a FunctionOperator for the similarity matrix $W$ that computes $Wx$, but in the paper we have $∇G = (W-I)^T(W-I)$.
I wanted to have a unified code where I can simply switch $∇G$ (L2 penalty, Huber loss, similarity regularization ...), but if you think QR factorization or a GMRES is the best way to do it, I will check that.

[ChrisRackauckas]

For that case, a truncated SVD can be made equivalent and will be more numerically stable: https://link.springer.com/article/10.1007/BF01937276 or matrix-free with GMRES via https://www.naturalspublishing.com/files/published/3v373jw5ql883r.pdf