Upgrade to v2.0

[andreasvarga]

[andreasvarga]

I updated my package MatrixEquations.jl to v2.0 and all tests involving LinearOperators failed! It is probably a trivial issue, but in this moment I would appreciate any hint on what I have to change when using the new version. A few days ago, I registered MatrixEquations.jl with the last version available before v2.0 and everything was OK.

[geoffroyleconte]

Hi! With v2.0 we updated the way to create operators, so that it is generic with matrix-vector products using mul! from LinearAlgebra. You can find an example here: https://juliasmoothoptimizers.github.io/LinearOperators.jl/stable/tutorial/#Using-functions . The prod! function should perform the operation: res = α * op * v + β * res.

We also removed PreallocatedLinearOperators, because you can now apply mul! to LinearOperators directly.

I hope that it helps, if it does not you can send me your operator and I'll try to help you update it to v2.0.

[andreasvarga]

Thanks. I wonder if this will be a simple exercise for me. Here is an example for the Lyapunov operator:

"""
    L = lyapop(A; disc = false, her = false) 

Define, for an `n x n` matrix `A`, the continuous Lyapunov operator `L:X -> AX+XA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-X` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is 
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated matrix `M = Matrix(L)` is 
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function lyapop(A; disc = false, her = false)
  n = LinearAlgebra.checksquare(A)
  T = eltype(A)
  function prod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A') - X)
      else
        Y = A * X
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        Y = A*X*A' - X
      else
        Y = A*X + X*A'
      end
      return Y[:]
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - X)
      else
        Y = X * A
        return triu2vec(Y + adjoint(Y))
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
         Y = adjoint(A)*X*A - X
       else
         Y = adjoint(A)*X + X*A
       end
       return Y[:]
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - X)
      else
        Y = X * A
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        return (A'*X*A - X )[:]
      else
        return (A'*X + X*A)[:]
      end
    end
  end
  her ? N = Int(n*(n+1)/2) : N = n*n
  return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
 

I have 16 operators defined and this was the simplest one! I wonder if there is an automatic safe way to perform the transitions. Is any compat tool available to suport the previous version? It would be great if you can keep also the previous way to define operators.

And, I must say frankly, I see no gain for my problems, which basically estimate the condition numbers of different kind of operators. The mul! function, if not carrefully implemented (e.g., to perform no additional operations for α = 1, β = 0), increases unnecessarily the number of operations, which is not desirable. So, the effort for users increases substantially (this is my understanding, please correct me if I am wrong).

I look forward to here your opinion.

[dpo]

@andreasvarga Keep in mind that there is no rush or necessity to upgrade at all. You can always restrict the version of LinearOperators required by your package in your Project.toml.

[geoffroyleconte]

I you want a quick fix you could try something like

function prod_op!(res, x, α, β)
  if β == 0
    res .= prod(x) .* α
  else
    res .= prod(x) .* α .+ β .* res
  end
end

with the prod function that you defined (and do something similar for the 2 other functions).

The main goal of v2 was to be able to call mul!(res, A, v, α , β) so that it is possible to write generic code when A is a matrix and A is an operator. To see performance gains, you should try to remove allocations from your prod function.

As you said, it increases the effort for the user, so we can try to add a feature that allows the user to keep his code the same in the future. For now, as @dpo said, if you do not want to make changes to your code, you can keep using v1.

[tknopp]

I find porting to v2 also quite difficult. The new opportunity to define 3-arg mult helps a lot here.

But in addition a found this helper function useful:

# Helper function to wrap a prod into a 5-args mul
function wrapProd(prod::Function)
  λ = (res, x, α, β) -> begin
    if β == zero(β)
      res .= prod(x) .* α
    else
      res .= prod(x) .* α .+ β .* res
    end
  end
  return λ
end

It can be used like this:

function WaveletOp(T::Type, shape, wt=wavelet(WT.db2))
  return LinearOperator{T}(prod(shape), prod(shape), false, false
            , wrapProd( x->vec( dwt(reshape(x,shape),wt) ) )
            , nothing
            , wrapProd( y->vec( idwt(reshape(y,shape),wt) ) ) )
end

I don't know if there are any performance issues but it makes porting much easier.