juliasmoothoptimizers / linearoperators.jl Goto Github PK

View Code? Open in Web Editor NEW

149.0 9.0 31.0 1.79 MB

Linear Operators for Julia

License: Other

Julia 100.00%

linear-operators linear-maps matrices preconditioner julia julia-language

linearoperators.jl's People

Contributors

Stargazers

Watchers

linearoperators.jl's Issues

Documentation outdated

Should be regenerated with modern tools.

License

On the documentation page

https://juliasmoothoptimizers.github.io/LinearOperators.jl/latest/

there is a GPLv3 badge but according to the license.md file this packages is under the MPL. Which is correct?

Block diagonal operator

It would be useful to be able to construct a block-diagonal operator from given blocks, e.g.,

P = blkdiag(M1, M2, ..., Mn)

would represent

[ M1           ]
[    M2        ]
[       ...    ]
[           Mn ]

Parametrized operators

Use a parametrized type instead of storing the data type inside the operator.

Unary operations on adjoint fail

I encoutered the following errors:

julia> L = LinearOperator(rand(2,2))
Linear operator
  nrow: 2
  ncol: 2
  eltype: Float64
  symmetric: false
  hermitian: false

julia> -L'
ERROR: type AdjointLinearOperator has no field nrow
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] -(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:260
 [3] top-level scope at REPL[56]:1

julia> 2*L'
ERROR: type AdjointLinearOperator has no field nrow
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] *(::Int64, ::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:311
 [3] top-level scope at REPL[57]:1

Implement `opHermitian` with preallocation

It's currently not possible to preallocate the output of an opHermitian. This could substantially speed up Krylov methods for symmetric and Hermitian systems.

Limited-memory quasi-Newton operators

limited-memory SR1 approximation
- in loop form, forward ~~and inverse~~: #14
- in compact form, forward ~~and inverse~~
- NB: inverse is trivial because SR1 is self dual
limited-memory BFGS approximation
- in compact form, forward and inverse
- damped update: #17 with corrections in #40
limited-memory of any approximation in the Broyden class
- in compact form, forward and inverse

Some operations involving adjoints are not yet supported

The following operations lead to errors:

julia> a = rand(2,2)
2×2 Array{Float64,2}:
 0.704653  0.99139 
 0.259309  0.292939

julia> A = LinearOperator(a)
Linear operator
  nrow: 2
  ncol: 2
  eltype: Float64
  symmetric: false
  hermitian: false

julia> A'+A
ERROR: type AdjointLinearOperator has no field symmetric
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] +(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:328
 [3] top-level scope at REPL[37]:1

julia> A'+a
ERROR: type AdjointLinearOperator has no field symmetric
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] +(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:328
 [3] +(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}, ::Array{Float64,2}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:334
 [4] top-level scope at REPL[46]:1

julia> [A A']
ERROR: type AdjointLinearOperator has no field nrow
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] hcat(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}, ::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:585
 [3] top-level scope at REPL[38]:1

julia> [A' A]
ERROR: type AdjointLinearOperator has no field nrow
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] hcat(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, 
Symbol("##12#15")){Array{Float64,2}}}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:585
 [3] top-level scope at REPL[39]:1

julia> [a  A']
ERROR: type AdjointLinearOperator has no field nrow
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] hcat(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}, ::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:585
 [3] hcat(::Array{Float64,2}, ::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:582
 [4] top-level scope at REPL[40]:1

julia> [A'; A]
ERROR: type AdjointLinearOperator has no field ncol
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] vcat(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, 
Symbol("##12#15")){Array{Float64,2}}}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:614
 [3] top-level scope at REPL[43]:1

julia> [A'; a]
ERROR: type AdjointLinearOperator has no field ncol
Stacktrace:
 [1] getproperty(::Any, ::Symbol) at .\Base.jl:20
 [2] vcat(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, 
Symbol("##12#15")){Array{Float64,2}}}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:614
 [3] vcat(::AdjointLinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, 
Symbol("##12#15")){Array{Float64,2}}}, ::Array{Float64,2}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:609
 [4] top-level scope at REPL[44]:1

Deprecation warning

Hi I get the following deprecation warning.

WARNING: importing deprecated binding Base.diag into LinearOperators.

LinearOperators concatenation functions have priority on matrices and fail

The following example has been executed after starting Julia:

julia> using LinearAlgebra

julia> a = rand(3,3)
3×3 Array{Float64,2}:
 0.909781  0.664048  0.36408
 0.408875  0.582388  0.863411
 0.392951  0.473908  0.463732

julia> [a' I]   # this works, the Base function hcat is used
3×6 Array{Float64,2}:
 0.909781  0.408875  0.392951  1.0  0.0  0.0
 0.664048  0.582388  0.473908  0.0  1.0  0.0
 0.36408   0.863411  0.463732  0.0  0.0  1.0

julia> using LinearOperators

julia> [a' I]   # this fails, the overloaded function hcat from LinearOperators is used
ERROR: MethodError: hcat(::Adjoint{Float64,Array{Float64,2}}, ::Array{Bool,2}) is ambiguous. Candidates:
  hcat(ops::Union{AbstractArray{T,2} where T, AbstractLinearOperator}...) in LinearOperators at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:601
  hcat(A::Union{AbstractArray{T,2}, AbstractArray{T,1}} where T...) in Base at abstractarray.jl:1313
Possible fix, define
  hcat(::Vararg{AbstractArray{T,2} where T,N} where N)
Stacktrace:
 [1] hcat(::Adjoint{Float64,Array{Float64,2}}, ::UniformScaling{Bool}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\uniformscaling.jl:320
 [2] top-level scope at REPL[6]:1

julia> [copy(a') I]   # this works again
3×6 Array{Float64,2}:
 0.909781  0.408875  0.392951  1.0  0.0  0.0
 0.664048  0.582388  0.473908  0.0  1.0  0.0
 0.36408   0.863411  0.463732  0.0  0.0  1.0

As it can be seen, the LinearOperators concatenation function hcat has priority on matrices and fails. It is any way to enforce the use of the function hcat from Base on matrices? Is this a Julia issue or a LinearOperators issue?

Better transposed operator

Currently, writing A' * u creates a new operator for A' so that preallocated storage for A cannot be taken advantage of. We need a proper transpose(A). This would also remove the need to use A.tprod(u) in certain Krylov.jl methods, which are not generic.

Recent changes

within the 0.5 something has changed which makes the following example not work anymore:

dft = LinearOperator{Float64}(10, 10, false, false,
                                     v -> fft(v),
                                     nothing,
                                     w -> ifft(w))

I get:

ERROR: MethodError: no method matching LinearOperator{Float64,F1,F2,F3} where F3<:Union{Nothing, Function} where F2<:Union{Nothing, Function} where F1<:Union{Nothing, Function}(::Int64, ::Int64, ::Bool, ::Bool, ::getfield(Main, Symbol("##31#33")), ::Nothing, ::getfield(Main, Symbol("##32#34")))
Stacktrace:

This example is still in the docs.

Document operators with preallocation

This is kind of crucial but currently not documented:

A = rand(m, n)
Ap = Vector{Float64}(n)  # buffer
Atq = Vector{Float64}(m)  # buffer
op = LinearOperator(m, n, false, false, p -> A_mul_B!(Ap,  A, p), q -> At_mul_B!(Atq, A, q))

The same can be done with functions:

op = LinearOperator(m, n, false, false, p -> Aprod!(Ap, p), q -> Atprod!(Atq, q))

where Aprod!() returns Ap and Atprod!() returns Atq.

Docs error

release on .travis should be 0.5. Docs didn't build.

julia 0.4 docstring compatibility issue?

I'm new to Julia. using LinearOperators gave this message:
WARNING: both Docile and Base export "@doc"; uses of it in module LinearOperators must be qualified

How to define A_mul_B! function for a LinearOperator?

I cannot find a way to do allocation-free multiplication of the linear operator by a vector. Is this feature implemented?

[discussion] Could we remove prod and define * for specific types?

Having prod, tprod, ctprod and their types makes it difficult for the compiler to guess the types of the outputs, which leads to instability. An alternative could be not assuming that each operator will have prod, and instead define subtypes and * for each one.

Example:

Remove generic *(A, v);
Implement subtype MatricialLinearOperator;
Define *(A, v) for this subtype.

We could then define everything in term of subtypes:

ScalingLinearOperator for scalar times operator;
MultiplicationLinearOperator for opA * opB;
hvcatLinearOperator;
AdjointLinearOperator (See #104);

One downside would be implementing a lot of new types and their interactions. For instance, adjoint of an operator would depend on the operator.

Related to #97 and #104

Missing methods for AdjointLinearOperator

using LinearOperators
A = LinearOperator(rand(10,10))
- A' * ones(10)
ERROR: type AdjointLinearOperator has no field nrow

Huge error log

Why is there a huge error log for this:

A = sprand(10, 10, 0.5)
v = zeros(3)
op = LinearOperator(A)
v[2] = op[2,3]

log:

ERROR: MethodError: Cannot `convert` an object of type LinearOperator{Float64,getfield(LinearOperators, Symbol("##36#39"))
{LinearOperator{Float64,getfield(LinearOperators, Symbol("##36#39"))
{LinearOperator{Int64,getfield(LinearOperators, Symbol("##155#157"))
{Array{Int64,1}},getfield(LinearOperators, Symbol("##156#158"))
{Int64,Array{Int64,1}},getfield(LinearOperators, Symbol("##156#158"))
{Int64,Array{Int64,1}}},LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){SparseMatrixCSC{Float64,Int64}},getfield(LinearOperators, 
Symbol("##11#14")){SparseMatrixCSC{Float64,Int64}},getfield(LinearOperators, 
Symbol("##12#15")){SparseMatrixCSC{Float64,Int64}}}},getfield(LinearOperators, 
Symbol("##37#40")){LinearOperator{Float64,getfield(LinearOperators, 
Symbol("##10#13")){SparseMatrixCSC{Float64,Int64}},getfield(LinearOperators, 
Symbol("##11#14")){SparseMatrixCSC{Float64,Int64}},getfield(LinearOperators, 
Symbol("##12#15...
[many lines later]
...1}}}}}}} to an object of type Float64

Exponentially long precompilation times with hcat and LinearOperators

I have observed that precompilation times increase exponentially with the number of LinearOperator instances in a call to hcat.

This behavior was observed in Julia 0.7, 1.0.3, 1.1 but NOT on Julia 0.6.4 on Mac OS 10.13.6 and on Julia 1.2 on Linux x64. I am using LinearOperators v0.5.4

Here is sample code:

# put in file lops.jl
# LinearOperator and hcat precompilation problem
using LinearOperators
N = 16; n = 5 
lops = [ LinearOperator(randn(n,n)) for i=1:N ]
C = hcat(lops...)
@time y = C*randn(n*length(lops))

On my machine include("lops.jl") takes 24 seconds on a first run, and then it drops to a more manageable 0.000036 seconds per run. The precompilation times I have observed by varying the size N of the hcat list (each timed on a fresh Julia 1.0 session) are:

N	time (seconds)
13	1.16
14	2.9
15	8.6
16	24
17	73
18	219

I should emphasize that the excessive overhead goes away on subsequent calls.

I think this is linked specifically to the hcat implementation in LinearOperators.jl since the following code achieves the same thing in a less elegant way but with a more reasonable 0.20 seconds time on first run (and also a trivial amount of time afterwards)

# put in file lops_fast.jl
# doing same operation as in lops.jl without hcat
using LinearOperators
N = 16; n = 5 
lops = [ LinearOperator(randn(n,n)) for i=1:N ]
function amul(x)
 y = fill(0.0,n)
 for i=1:N
  y += lops[i]*x[n*(i-1) .+ (1:n)]
 end
 return y
end
function amult(y)
 x = fill(0.0,n*N)
 for i=1:N
  x[n*(i-1) .+ (1:n)] = lops[i]'*y
 end
 return x
end
C = LinearOperator(Float64,n,N*n,false,false,amul,amult,amult)
@time y = C*randn(n*length(lops))

Refactor source code

Split source files into something like "utilities.jl", "types.jl", "basics.jl", "operations.jl", "specials.jl".

Original: #56

Indexing

Given an operator A, it would be useful to be able to index it, e.g., A[1:2:9, cols]. An elegant way of doing this could be R * A * E where R is a RestrictionOperator and E is an ExtensionOperator.

norm(matrix) -> opnorm(matrix) in Julia 0.7

It looks like may be using norm(matrix). In Julia 0.7, this will compute the Frobenius norm (vecnorm in Julia 0.6), due to JuliaLang/julia#27401. If you want the induced/operator norm as in Julia 0.6, use opnorm(matrix) instead, or Compat.opnorm(matrix) to work in 0.6 and 0.7 (JuliaLang/Compat.jl#577).

Note that, for testing purposes, rather than @test norm(A - B) ≤ tol, it is usually preferred to do @test A ≈ B or @test A ≈ B rtol=... (which uses isapprox).

Commutation with numbers and operators

* between operators and Numbers is assumed to be a commutative operation, but this could not be the expected behaviour for some specific Numbers (quaternions? ).

check_positive_definite semi kwarg not working

In LinearOperators.jl, line 496, I had to change the line to
check_positive_definite(M :: AbstractMatrix; semi=false) = check_positive_definite(LinearOperator(M); semi=false)

in order to use the semi keyword argument

Maybe something to check

Some operations do not work although they (probably) should

The following code illustrates this issue:

M = ones(Rational,3,3)
v = ones(Int,3)
M*v  # OK, the result is a rational vector
LinearOperator(M)*v # failure

The execution of this code produces:


ERROR: TypeError: in typeassert, expected Array{Real,1}, got a value of type Array{Rational,1}
Stacktrace:
 [1] *(::LinearOperator{Rational}, ::Array{Int64,1}) at C:\Users\Andreas\.julia\packages\LinearOperators\z6AM5\src\operations.jl:7
 [2] top-level scope at REPL[126]:1
 [3] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088

I wonder what is the right approach in determining the type of the elements of the result of an operation M*v. For example, in my concrete case, M is an inverse operator to solve a matrix equation, with the underlying matrices of a certain element type T (say T = Float32). The vector v can have any element type T1 (say T1 = Float64) such that it can be still converted to the element type T, without error. In this case, the result M*v would have element type T, which would be, for me, the desirable type. Still note that for M a matrix, the result of M*v has element type T (and not T1), which involves a conversion of M (generally not desirable in my case).

Some (many?) proposals

Hello, newbie here. I have in mind some features and possible optimization that
could help to improve this package, also in terms of maintainability. I wanted
to open this issue for having a feedback from the maintainers before starting
to open Pull Requests.

My ideas consist of

Splitting source files into something like "utilities.jl", "types.jl",
"basics.jl", "operations.jl", "specials.jl". I think that each file name
is self-explicative.
Making functions for AbstractLinearOperator more generic, substituting
explicit fields invocations (e.g. op.symmetric) with traits (e.g.
issymmetric(op)).
Making LinearOperator immutable (don't have find any function nor
reason why it should be mutable but, please, correct me if I am wrong).
Parametrizing LinearOperator{T} with the type of transformation
function, (to become LinearOperator{T, Fprod, Ftprod, Fctprod}) so that
there is more specialization and in theory a better performance on code
(this has to be tested).
Using undef to signal non-specified tprods and ctprods instead
of nothing, because in my opinion it would be more self-describing.
Creating a new struct to replace v -> M * v, something like a callable
```
struct MatrixProduct{T <: AbstractMatrix}
    M::T
end
(mp::MatrixProduct)(v) = mp.M * v
```
to preserve information about M and use it to combine linear operators in
a more efficient manner.
Removing information about being symmetric and being hermitian from
fields and creating two new AbstractLinearOperators such as
SymmetricLinearOperator and HermitianLinearOperator with which to
dispatch, instead of using many if-else statements. Symmetry can then
be checked as a trait and it could be useful a function linearoperator(x)
that transform a Matrix or a function in the most informative type.

I am ready to take care of these PRs. None of them seems breaking to me,
but first of all let me know what do you think about them.

In-place solve for inverse operators

I would be nice to add in-place solves to LinearOperators in the linalg.jl file.

PreallocatedLinearOperators with Int types

A = rand(Int, 10, 10)
b = rand(10)
M = PreallocatedLinearOperator(A)

julia> M * b
ERROR: InexactError: Int64(-9.58046120180543e18)
Stacktrace:
 [1] Type at ./float.jl:703 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:767 [inlined]
 [4] generic_matvecmul!(::Array{Int64,1}, ::Char, ::Array{Int64,2}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:549
 [5] mul! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:76 [inlined]
 [6] #2 at /home/alexis/.julia/dev/LinearOperators/src/PreallocatedLinearOperators.jl:62 [inlined]
 [7] *(::PreallocatedLinearOperator{Int64,getfield(LinearOperators, Symbol("##2#5")){Array{Int64,1},Array{Int64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Int64,1},Array{Int64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Int64,1},Array{Int64,2}}}, ::Array{Float64,1}) at /home/alexis/.julia/dev/LinearOperators/src/LinearOperators.jl:227
 [8] top-level scope at none:0

Maybe we should preallocated Int matrix with Float vectors?
M * b will only work if b <: AbstractVector{Integer}, it's a bit annoying.
I think we can do the following modifiication:

function PreallocatedLinearOperator(M :: AbstractMatrix{T};
                                    symmetric=false, hermitian=false) where T
  nrow, ncol = size(M)
  local Mv, Mtu, Maw
   if T <: Integer
    if symmetric
      Maw = Mtu = Mv = Vector{Float64}(undef, nrow)
    else
      Mv = Vector{Float64}(undef, nrow)
      Maw = Mtu = Vector{Float64}(undef, ncol)
    end
  elseif T <: AbstractFloat
    if symmetric
      Maw = Mtu = Mv = Vector{T}(undef, nrow)
    else
      Mv = Vector{T}(undef, nrow)
      Maw = Mtu = Vector{T}(undef, ncol)
    end
  else
    if symmetric && hermitian # Actually real, but T is not
      Maw = Mtu = Mv = Vector{T}(undef, nrow)
    elseif symmetric
      Mtu = Mv = Vector{T}(undef, nrow)
      Maw = Vector{T}(undef, ncol)
    elseif symmetric
      Mv = Vector{T}(undef, nrow)
      Maw = Mtu = Vector{T}(undef, ncol)
    else
      Mv = Vector{T}(undef, nrow)
      Mtu = Vector{T}(undef, ncol)
      Maw = Vector{T}(undef, ncol)
    end
  end
  PreallocatedLinearOperator(Mv, Mtu, Maw, M, symmetric=symmetric, hermitian=hermitian)
end

Transpose syntax deprecated in 1.0

The transpose syntax appears to be broken in 1.0. Here is a minimum working example:

using LinearOperators; A = LinearOperator(rand(10,10)); A.'

This produces a syntax error:
ERROR: syntax: the ".'" operator is discontinued

$ julia --version
julia version 1.0.1

real kron operator applied to complex vector fails

The following sequence leads to an error:

julia> a = rand(2,2)
2×2 Array{Float64,2}:
 0.721074  0.441174
 0.902213  0.442013

julia> A = LinearOperator(a)
Linear operator
  nrow: 2
  ncol: 2
  eltype: Float64
  symmetric: false
  hermitian: false

julia> x = rand(Complex{Float64},4)
4-element Array{Complex{Float64},1}:
 0.12261904045152261 + 0.7466118601168086im
  0.3218599517061915 + 0.7258818549661219im
   0.701833424057472 + 0.2956429356370669im
 0.47354227893318335 + 0.31005578035910797im

julia> kron(a,A)*x
ERROR: InexactError: Float64(0.12261904045152261 + 0.7466118601168086im)
Stacktrace:
 [1] Type at .\complex.jl:37 [inlined]
 [2] convert at .\number.jl:7 [inlined]
 [3] setindex! at .\array.jl:766 [inlined]
 [4] copyto!(::IndexLinear, ::Array{Float64,1}, ::IndexLinear, ::Array{Complex{Float64},1}) at .\abstractarray.jl:807
 [5] copyto! at .\abstractarray.jl:799 [inlined]
 [6] Type at .\array.jl:482 [inlined]
 [7] convert(::Type{Array{Float64,1}}, ::Array{Complex{Float64},1}) at .\array.jl:474
 [8] (::getfield(LinearOperators, Symbol("#prod#139")){LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},Int64,Int64,DataType})(::Array{Complex{Float64},1}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\kron.jl:15
 [9] *(::LinearOperator{Float64,getfield(LinearOperators, Symbol("#prod#139")){LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},Int64,Int64,DataType},getfield(LinearOperators, Symbol("#tprod#140")){LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},Int64,Int64,DataType},getfield(LinearOperators, Symbol("#ctprod#141")){LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},LinearOperator{Float64,getfield(LinearOperators, Symbol("##10#13")){Array{Float64,2}},getfield(LinearOperators, Symbol("##11#14")){Array{Float64,2}},getfield(LinearOperators, Symbol("##12#15")){Array{Float64,2}}},Int64,Int64,DataType}}, ::Array{Complex{Float64},1}) at C:\Users\Andreas\.julia\packages\LinearOperators\La2Ey\src\LinearOperators.jl:229
 [10] top-level scope at REPL[50]:1

In contrast, the following works:

julia> kron(a,a)*x
4-element Array{Complex{Float64},1}:
  0.481579783754581 + 0.7735129831868585im
 0.5540519254548673 + 0.8952113798255193im
 0.5239164901392808 + 0.9293333332436311im
 0.6005678082785425 + 1.0756846076728297im

Fixing the error is simple, for example by using the sequence:

  function prod(x)
     T1 = promote_type(T, eltype(x))
     X = reshape(convert(Vector{T1}, x), q, n)
    return Matrix(B * X * transpose(A))[:]
  end

Implement opKron

Can't use the full function

I'm not sure if I haven't searched correctly but I can't use the full function to materialize function while calling Matrix(op) does work. Since there is no implementation that I could find for full I reckon its just a documentation issue.

Warning messages

Why do we see warning messages of the form

WARNING: deprecated syntax "function .*(...)".
Use "function Base.broadcast(::typeof(*), ...)" instead.

despite the fact that .* is only declared for Julia prior to 0.6?

General preallocated operators

Currently, preallocated linear operators are only defined for matrix input.

New release

After #29 ?

TypeError

julia> B = LBFGSOperator(5);

julia> x = rand(10);

julia> B * view(x, 1:5)
ERROR: TypeError: in typeassert, expected SubArray{Float64,1,P,I,L} where L where I where P, got a value of type Array{Float64,1}
Stacktrace:
 [1] *(::LBFGSOperator{Float64}, ::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}) at /Users/dpo/dev/JSO/LinearOperators.jl/src/operations.jl:7
 [2] top-level scope at REPL[54]:1

PreallocatedOperators (name pending)

cc @dpo @Xx-fractale-xX

For each current operator opExample:

create opExample!(W, ...) which uses W as workspace for opExample operation;
if possible, change opExample(...) to create workspace W and use above opExample!(W, ...).

Matrix-vector products with CuArrays and SparseArrays

julia> A = cu(rand(3, 3))
3×3 CuArray{Float32,2,Nothing}:
 0.678375  0.0608932  0.571121
 0.639773  0.5497     0.642157
 0.432494  0.744255   0.839123

julia> b = cu(rand(3))
3-element CuArray{Float32,1,Nothing}:
 0.42534474 
 0.119808994
 0.5397634  

julia> op = LinearOperator(A)
Linear operator
  nrow: 3
  ncol: 3
  eltype: Float32
  symmetric: false
  hermitian: false
  nprod:   0
  ntprod:  0
  nctprod: 0

julia> op * b
ERROR: TypeError: in typeassert, expected Array{Float32,1}, got CuArray{Float32,1,Nothing}
Stacktrace:
 [1] *(::LinearOperator{Float32}, ::CuArray{Float32,1,Nothing}) at /home/montalex/.julia/packages/LinearOperators/WaECw/src/operations.jl:7
 [2] top-level scope at REPL[6]:1

Add counters for number of products performed

Linear operators acting on symmetric matrices

I implemented several linear operators related to solve control problems. An example is the Lyapunov operator L: X -> AX+XA', where A is a squre matrix. This operator acts usually on symmetric/hermitian matrices X and is a special case of the more general Sylvester operator S: X -> AX+XB, with A and B square matrices (not necessarily of the same size). The definition of L as a linear operator is possible regardless X is symmetric/hermitian or not and Matrix(L) is the usual Kronecker-expansion based matrix. However, the definition of the inverse operator inv(L) : Y -> X involves the solution of a Lyapunov equation AX+XA' + Y = 0, where for a symmetric/hermitian Y the resulting X is also symmetric/hermitian. The solvers of Lyapunov equations usually exploit the symmetry of the solutions, by computing, for example, only the upper triangular part of X (the lower triangular part results via symmetry). It is possible to define inv(L) ignoring the symmetry in the input data. Unfortunately, in this case, some functions in the LinearOperators collection will fail, as for example Matrix(inv(L)), which will not generate the inverse operator matrix due to the restriction on the symmetry of the specialized solvers. To cope with this aspect, I was forced to use the more general solvers for Sylvester equations to compute the full solution X, with the associated efficiency losses.

I wonder if the problem of restrining the domain of input data to form the products as L * vec(X) can be addressed somehow, by assuming certain structural constraints on X (e.g., symetric, or hermitian or even diagonal).

Many thanks in advance for your time considering this question.

Note: The implemented linear operators belong to the recently developed (not yet registered) package MatrixEquation.jl.

Can't multiply with views

Operator-vector product breaks down if the vector is a view:

using LinearOperators

A = ones(2, 2)
opA = LinearOperator(A)

x = ones(2)
@views xv = x[:]

opA * x   # works fine
opA * xv  # error

raises the error

ERROR: TypeError: in typeassert, expected SubArray{Float64,1,P,I,L} where L where I where P, got a value of type Array{Float64,1}
Stacktrace:
 [1] *(::LinearOperator{Float64}, ::SubArray{Float64,1,Array{Float64,1},Tuple{Base.Slice{Base.OneTo{Int64}}},true}) at /home/mtanneau/.julia/packages/LinearOperators/z6AM5/src/operations.jl:7
 [2] top-level scope at REPL[7]:1

Make Linear Operators immutable

Original: #56

New ambiguity on Julia 1.3

While running the package tests for the upcoming 1.3, a new ambiguity in "HybridSystems.jl" was detected but it seems to originate from this package:

Cholesky and LDL: Error During Test at /root/.julia/packages/LinearOperators/M4IRn/test/test_linop.jl:369
  Got exception outside of a @test
  MethodError: hvcat(::Tuple{Int64,Int64}, ::Array{Float64,2}, ::Adjoint{Float64,Array{Float64,2}}, ::Array{Float64,2}, ::Array{Float64,2}) is ambiguous. Candidates:
    hvcat(rows::Tuple{Vararg{Int64,N} where N}, xs::Union{AbstractArray{T,1}, AbstractArray{T,2}}...) where T in Base at abstractarray.jl:1634
    hvcat(rows::Tuple{Vararg{Int64,N} where N}, ops::Union{AbstractArray{T,2} where T, AbstractLinearOperator}...) in LinearOperators at /root/.julia/packages/LinearOperators/M4IRn/src/LinearOperators.jl:697
    hvcat(rows::Tuple{Vararg{Int64,N} where N}, xs::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T...) in Base at abstractarray.jl:1633
  Possible fix, define
    hvcat(::Tuple{Vararg{Int64,N} where N}, ::Vararg{Union{AbstractArray{T,1}, AbstractArray{T,2}},N} where N)

We believe this is a real ambiguity that got detected due to improvements in the Julia type system and thus needs to be fixed in this package.

Matrices in block operators

julia> K = [ opEye(10) opZeros(10,5) ; opZeros(5,10) speye(5) ]
ERROR: DimensionMismatch("mismatch in dimension 1 (expected 1 got 5)")
Stacktrace:
 [1] _cs(::Int64, ::Bool, ::Int64, ::Int64) at ./abstractarray.jl:1193
 [2] _cshp at ./abstractarray.jl:1189 [inlined]
 [3] cat_shape(::Tuple{Bool,Bool}, ::Tuple{Int64,Int64}, ::Tuple{Int64,Int64}) at ./abstractarray.jl:1170 (repeats 2 times)
 [4] cat_t(::Type{T} where T, ::Type{T} where T, ::LinearOperators.LinearOperator{Float64}, ::Vararg{Any,N} where N) at ./abstractarray.jl:1203
 [5] hcat(::LinearOperators.LinearOperator{Float64}, ::SparseMatrixCSC{Float64,Int64}) at ./abstractarray.jl:1301
 [6] hvcat(::Tuple{Int64,Int64}, ::LinearOperators.LinearOperator{Float64}, ::Vararg{Any,N} where N) at ./abstractarray.jl:1477

Matrix(op) may generate incorrect result if op * x changes x

Matrix(op) may generate incorrect result in the case when op * x changes the values stored in x. Typical example is when op is a linear equation solver which overwrites the result in the provided right hand side. Then the repeated computation of op * ei in the following code may lead to a change of ei after each execution of op * ei. A solution would be to move

ei = zeros(eltype(op), n)

inside the loop. Otherwise, the change of x when computing op * x must be explicitly forbiden.

function Base.Matrix(op :: AbstractLinearOperator)
  (m, n) = size(op)
  A = Array{eltype(op)}(undef, m, n)
  ei = zeros(eltype(op), n)
  for i = 1 : n
    ei[i] = 1
    A[:, i] = op * ei
    ei[i] = 0
  end
  return A
end

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Missing methods

julia> using LinearAlgebra, LinearOperators

julia> A = randn(500,500); op = LinearOperator(A);

julia> y = zeros(500); x = ones(500);

julia> mul!(y, op, x)
ERROR: MethodError: no method matching length(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}})
Closest candidates are:
  length(::Core.SimpleVector) at essentials.jl:582
  length(::Base.MethodList) at reflection.jl:732
  length(::Core.MethodTable) at reflection.jl:806
  ...
Stacktrace:
 [1] _similar_for(::UnitRange{Int64}, ::Type, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}, ::Base.HasLength) at ./array.jl:532
 [2] _collect(::UnitRange{Int64}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}, ::Base.HasEltype, ::Base.HasLength) at ./array.jl:563
 [3] collect(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}) at ./array.jl:557
 [4] broadcastable(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}) at ./broadcast.jl:617
 [5] broadcasted at ./broadcast.jl:1171 [inlined]
 [6] mul!(::Array{Float64,1}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}, ::Array{Float64,1}) at /Users/dpo/.julia/packages/LinearOperators/aODV2/src/LinearOperators.jl:295
 [7] top-level scope at none:0

julia> import Base.length

julia> length(op::AbstractLinearOperator) = op.nrow * op.ncol
length (generic function with 83 methods)

julia> mul!(y, op, x)
ERROR: MethodError: no method matching iterate(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}})
Closest candidates are:
  iterate(::Core.SimpleVector) at essentials.jl:589
  iterate(::Core.SimpleVector, ::Any) at essentials.jl:589
  iterate(::ExponentialBackOff) at error.jl:171
  ...
Stacktrace:
 [1] copyto!(::Array{Float64,1}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}) at ./abstractarray.jl:646
 [2] _collect(::UnitRange{Int64}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}, ::Base.HasEltype, ::Base.HasLength) at ./array.jl:563
 [3] collect(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}) at ./array.jl:557
 [4] broadcastable(::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}) at ./broadcast.jl:617
 [5] broadcasted at ./broadcast.jl:1171 [inlined]
 [6] mul!(::Array{Float64,1}, ::LinearOperator{Float64,getfield(LinearOperators, Symbol("##2#5")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##3#6")){Array{Float64,1},Array{Float64,2}},getfield(LinearOperators, Symbol("##4#7")){Array{Float64,1},Array{Float64,2}}}, ::Array{Float64,1}) at /Users/dpo/.julia/packages/LinearOperators/aODV2/src/LinearOperators.jl:295
 [7] top-level scope at none:0

Compatibility with Julia 0.4+

Nothing has become Void in 0.4.

juliasmoothoptimizers / linearoperators.jl Goto Github PK

linearoperators.jl's People

Contributors

Stargazers

Watchers

Forkers

linearoperators.jl's Issues

Recommend Projects

Recommend Topics

Recommend Org