Inverse problems

[jagot]

I managed to get some of the basis inversions to work, but not all:

using FiniteDifferencesQuasi
using LinearAlgebra
using Test

R = FiniteDifferences(20,0.2)
R⁻¹ = pinv(R)
# # Actually, we want
# R⁻¹ = inv(R)
# # but then we get
# # ERROR: LoadError: MethodError: no method matching ContinuumArrays.AlephInfinity{1}(::Int64)

u = R*rand(size(R,2))
v = R*rand(size(R,2))

uvc = [u.args[2] v.args[2]]
uv = R*uvc

ut = u'

@test R⁻¹*R == I
@test R*R⁻¹ == I
@test R⁻¹*u === u.args[2]
@test R⁻¹*uv === uvc

# # We also want this to work
# @test ut*R⁻¹' === ut.args[1]
# # ERROR: MethodError: no method matching axes(::UniformScaling{Float64}, ::Int64)

M = Diagonal(2ones(size(R,2)))
# # Finally, we'd like to write operators as
# O = R*M*R⁻¹
# # but alas
# # ERROR: LoadError: MethodError: no method matching ContinuumArrays.AlephInfinity{1}(::Int64)
# @test R⁻¹*O*u == M*R⁻¹*u

I've not yet figured out why inv does not work, to me it seems most things are mirrored compared to pinv in e.g. inv.jl.

[dlfivefifty]

Please paste the whole error message

[jagot]

Sorry about that:

julia> inv(R)
ERROR: MethodError: no method matching ContinuumArrays.AlephInfinity{1}(::Int64)
Closest candidates are:
  ContinuumArrays.AlephInfinity{1}(::T<:Number) where T<:Number at boot.jl:718
  ContinuumArrays.AlephInfinity{1}() where N at /Users/jagot/.julia/packages/ContinuumArrays/8Way3/src/ContinuumArrays.jl:25
  ContinuumArrays.AlephInfinity{1}(::Float16) where T<:Integer at float.jl:71
  ...
Stacktrace:
 [1] convert(::Type{ContinuumArrays.AlephInfinity{1}}, ::Int64) at ./number.jl:7
 [2] zero(::Type{ContinuumArrays.AlephInfinity{1}}) at ./number.jl:266
 [3] Base.OneTo{ContinuumArrays.AlephInfinity{1}}(::ContinuumArrays.AlephInfinity{1}) at ./range.jl:309
 [4] Base.OneTo(::ContinuumArrays.AlephInfinity{1}) at ./range.jl:318
 [5] axes1(::QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}) at ./abstractarray.jl:96
 [6] map(::typeof(Base.axes1), ::Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}}) at ./tuple.jl:140
 [7] axes(::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndices{2,Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}},Tuple{Int64,Float64}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/multidimensional.jl:308
 [8] _array_for(::Type{Float64}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndices{2,Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}},Tuple{Int64,Float64}}, ::Base.HasShape{2}) at ./array.jl:598
 [9] Array{Float64,N} where N(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/abstractquasiarray.jl:20
 [10] Array(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/abstractquasiarray.jl:22
 [11] QuasiArrays.QuasiArray(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/quasiarray.jl:35
 [12] QuasiArrays.QuasiArray(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/lazyquasiarrays.jl:82
 [13] copy(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/lazyquasiarrays.jl:81
 [14] materialize(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(inv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/lazyapplying.jl:40
 [15] apply at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/lazyapplying.jl:37 [inlined]
 [16] inv(::FiniteDifferences{Float64,Int64}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:59
 [17] top-level scope at REPL[9]:1

julia> ut*R⁻¹'
ERROR: MethodError: no method matching axes(::UniformScaling{Float64}, ::Int64)
Closest candidates are:
  axes(::Core.SimpleVector, ::Integer) at essentials.jl:609
  axes(::Number, ::Integer) at number.jl:65
  axes(::ND<:FiniteDifferencesQuasi.NumerovDerivative, ::Any...) where ND<:FiniteDifferencesQuasi.NumerovDerivative at /Users/jagot/.julia/dev/FiniteDifferencesQuasi/src/FiniteDifferencesQuasi.jl:377
  ...
Stacktrace:
 [1] _check_mul_axes(::Adjoint{Float64,Array{Float64,1}}, ::UniformScaling{Float64}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/linalg/mul.jl:19
 [2] check_mul_axes(::Adjoint{Float64,Array{Float64,1}}, ::UniformScaling{Float64}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/linalg/mul.jl:21
 [3] check_applied_axes(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},UniformScaling{Float64}}}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/linalg/mul.jl:25
 [4] instantiate(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},UniformScaling{Float64}}}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/lazyapplying.jl:22
 [5] materialize(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},UniformScaling{Float64}}}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/lazyapplying.jl:40
 [6] fullmaterialize(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},UniformScaling{Float64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:99
 [7] fullmaterialize(::LazyArrays.Applied{QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},QuasiArrays.QuasiAdjoint{Float64,FiniteDifferences{Float64,Int64}},QuasiArrays.QuasiAdjoint{Float64,QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:114
 [8] fullmaterialize at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:120 [inlined]
 [9] *(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},QuasiArrays.QuasiAdjoint{Float64,FiniteDifferences{Float64,Int64}}}}, ::QuasiArrays.QuasiAdjoint{Float64,QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:55
 [10] top-level scope at REPL[16]:1

julia> R*M*R⁻¹
ERROR: MethodError: no method matching ContinuumArrays.AlephInfinity{1}(::Int64)
Closest candidates are:
  ContinuumArrays.AlephInfinity{1}(::T<:Number) where T<:Number at boot.jl:718
  ContinuumArrays.AlephInfinity{1}() where N at /Users/jagot/.julia/packages/ContinuumArrays/8Way3/src/ContinuumArrays.jl:25
  ContinuumArrays.AlephInfinity{1}(::Float16) where T<:Integer at float.jl:71
  ...
Stacktrace:
 [1] convert(::Type{ContinuumArrays.AlephInfinity{1}}, ::Int64) at ./number.jl:7
 [2] zero(::Type{ContinuumArrays.AlephInfinity{1}}) at ./number.jl:266
 [3] Base.OneTo{ContinuumArrays.AlephInfinity{1}}(::ContinuumArrays.AlephInfinity{1}) at ./range.jl:309
 [4] Base.OneTo(::ContinuumArrays.AlephInfinity{1}) at ./range.jl:318
 [5] axes1(::QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}) at ./abstractarray.jl:96
 [6] map(::typeof(Base.axes1), ::Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}}) at ./tuple.jl:140
 [7] axes(::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndices{2,Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}},Tuple{Int64,Float64}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/multidimensional.jl:308
 [8] _array_for(::Type{Float64}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndices{2,Tuple{Base.OneTo{Int64},QuasiArrays.Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}}},Tuple{Int64,Float64}}, ::Base.HasShape{2}) at ./array.jl:598
 [9] Array{Float64,N} where N(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/abstractquasiarray.jl:20
 [10] Array(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/abstractquasiarray.jl:22
 [11] QuasiArrays.QuasiArray(::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/quasiarray.jl:35
 [12] QuasiArrays.QuasiArray(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/lazyquasiarrays.jl:82
 [13] copy(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/lazyquasiarrays.jl:81
 [14] materialize(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/LazyArrays/znbQW/src/lazyapplying.jl:40
 [15] fullmaterialize(::LazyArrays.Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(*),Tuple{Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:99
 [16] fullmaterialize(::LazyArrays.Applied{QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{FiniteDifferences{Float64,Int64},Diagonal{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:113
 [17] fullmaterialize at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:120 [inlined]
 [18] *(::FiniteDifferences{Float64,Int64}, ::Diagonal{Float64,Array{Float64,1}}, ::QuasiArrays.ApplyQuasiArray{Float64,2,typeof(pinv),Tuple{FiniteDifferences{Float64,Int64}}}) at /Users/jagot/.julia/packages/QuasiArrays/3525V/src/matmul.jl:54
 [19] top-level scope at REPL[17]:1

[dlfivefifty]

Ok we need ApplyStyle(pinv, typeof(R)) to return LazyQuasiArrayApplyStyle() to tell it not to try to materialize. Will check why this isn't already the case.

[jagot]

Gotcha. The same error occurs when trying form the (pseudo-)inverse of a restricted basis: inv(R[:,2:end-1]).