jp/deflation

[ioannisPApapadopoulos]

Here are the necessary modifications to the Newton iteration in the equilibriummeasure function to utilize deflation when looking for multiple solutions of a one-interval log-kernel equilbrium measures problem. The Newton update is computed as normal and then (nonlinearly) scaled to prevent the Newton method converging to an already discovered solution.

An example of a problem with 9 solutions (3 feasible) is found in examples/non-unique-equilibrium-measure.jl.

[TSGut]

This is really cool!

In terms of the terminology we've been using, are they all equilibrium measures? I think the global minimizer is called the equilibrium measure and the other potential local minimizers are just referred to as steady states and non uniqueness of steady states is not necessarily the same as non uniqueness of the equilibrium measures. We should maybe compute the energies along with the measures to check this.

[ioannisPApapadopoulos]

great, thanks! in terms of the code, does that mean just computing the energies of all the 9 solutions in the example script? Is there anything else you had in mind?

[dlfivefifty]

Yes, the energies should be greater than the equilibrium measure which has 3 intervals, and can currently only be computed in Mathematica

you’ll need the log transform implemented to compute the energy, though point wise evaluation is fine and is easier to compute. Check out StieltjesPoint in stieltjes.jl for an example with the Cauchy/Stieltjes kernel

[ioannisPApapadopoulos]

Am I wrong in thinking that implementing the log transform for weighted U in this example can be done by simply re-implementing the already coded up log transform for weighted T (stieltjes.jl, lines 157-171)?

[dlfivefifty]

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/94d58435970c692060351b68a81d08b41322ff74/src/stieltjes.jl#L157

is only for evaluating the log transform on the same interval. You can tell because LogKernel{T,R,D} where R is the type of the range interval and D is the type of the domain interval, in this case R == D == ChebyshevInterval, that is we know the domain and range are the same.

What we need for two-intervals is evaluating when R is a different interval.

For comparison, for the Hilbert/Stieltjes case we have R == D == ChebyshevInterval here:

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/94d58435970c692060351b68a81d08b41322ff74/src/stieltjes.jl#L87

and R ≠ D == ChebyshevInterval here:

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/94d58435970c692060351b68a81d08b41322ff74/src/stieltjes.jl#L256

[ioannisPApapadopoulos]

ok! and here we only want to perform the log transform over the range where the (feasible) measures are supported? So not on the whole domain (a,b)?

[dlfivefifty]

We want two things:

1. Evaluation of the energy at a point. That is, we want to construct the functional/row vector corresponding to

x = 0.1 # some real point, either on the domain or not
W = Weighted(ChebyshevU())
t = axes(W,1)
log.(abs.(x .- t')) * W # should return a row-vector

This could be done by creating

const LogKernelPoint{T,V,D} = BroadcastQuasiMatrix{T,typeof(log),BroadcastQuasiMatrix{T,typeof(abs),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{T,QuasiAdjoint{V,Inclusion{V,D}}}}}}}

by mimicking StieltjesPoint, see
https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/94d58435970c692060351b68a81d08b41322ff74/src/stieltjes.jl#L214
Note in Julia we want it to be a transpose or adjoint of a vector so that multiplication with coefficients gives a constant.

2. Re-expansion on another interval in chebyshevt(a..b). This isn't actually needed for now but is useful for solving other SIEs with log kernels on two intervals.

[ioannisPApapadopoulos]

So should the following already be working?

x = 0.1 # some real point, either on the domain or not
W = Weighted(ChebyshevU())
t = axes(W,1)
inv.(x .- t')' * W # should return a row-vector

I currently get a DimensionMismatch error.

DimensionMismatch("Second axis of A, Base.OneTo(1), and first axis of B, Inclusion(-1.0..1.0 (Chebyshev)) must match")

[dlfivefifty]

You have an extra ‘

[dlfivefifty]

because I did.. fixed my example

[ioannisPApapadopoulos]

without the extra ', I am still getting an error:

MethodError: no method matching materialize!(::ArrayLayouts.MulAdd{ArrayLayouts.UnknownLayout, ArrayLayouts.UnknownLayout, LazyArrays.PaddedLayout{ArrayLayouts.DenseColumnMajor}, Float64, QuasiArrays.SubQuasiArray{Float64, 2, QuasiArrays.PermutedDimsQuasiArray{Float64, 2, (2, 1), (2, 1), Weighted{Float64, ChebyshevU{Float64}}}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, false}, QuasiArrays.SubQuasiArray{Float64, 1, QuasiArrays.PermutedDimsQuasiArray{Float64, 2, (2, 1), (2, 1), QuasiArrays.BroadcastQuasiMatrix{Float64, typeof(inv), Tuple{QuasiArrays.BroadcastQuasiMatrix{Float64, typeof(-), Tuple{QuasiArrays.QuasiAdjoint{Float64, Inclusion{Float64, ChebyshevInterval{Float64}}}, Inclusion{Float64, ChebyshevInterval{Float64}}}}}}}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, Float64}, false}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}})
Closest candidates are:
  materialize!(::Any, ::Base.Broadcast.Broadcasted{BS, Axes, F, Args} where {Axes, F, Args<:Tuple}) where {Style, BS<:BlockArrays.AbstractBlockStyle} at C:\Users\john.papad\.julia\packages\BlockArrays\rGDmR\src\blockbroadcast.jl:147
  materialize!(::Any, ::Base.Broadcast.Broadcasted{Style, Axes, F, Args} where {Axes, F, Args<:Tuple}) where Style at broadcast.jl:890
  materialize!(::Any, ::Any) at broadcast.jl:886
  ...

[dlfivefifty]

Ah, for some reason it's only working for complex numbers right now:

julia> inv.(x .- t') * W
transpose((ComplexF64) .* ((ComplexF64) .^ ℵ₀-element InfiniteArrays.OneToInf{Int64} with indices OneToInf() with indices OneToInf()) with indices OneToInf()) with indices Base.OneTo(1)×OneToInf():
 0.0917376-1.29575im  -0.531751-0.0756742im  -0.0467394+0.21711im  …  

I'll have to look into that...

[dlfivefifty]

Actually the issue isn't real or complex, it's that the code for evaluation when x in ChebyshevInterval() is broken.

JuliaApproximation/ClassicalOrthogonalPolynomials.jl#66

has a fix.

[ioannisPApapadopoulos]

Thanks @dlfivefifty! Switching to the branch dl/offhilbertcolsupportbug did the trick and it works on my end now too!

I suppose that implementing the evaluation of the energy at a point for the log-kernel would be part of a seperate PR for ClassicalOrthogonalPolynomials.jl. Would you want to merge this deflation PR into EquilibriumMeasures.jl or leave it as a branch for now?

[dlfivefifty]

Let’s merge this, but perhaps it would be helpful to make it support other types like BigFloat . Do you want to have a try, or maybe at this point it is more helpful if I go through the code and make the necessary changes so you can learn what to do?

[ioannisPApapadopoulos]

I am willing to give it a go! This would be my first time dealing with BigFloat, do you have a pointer to somewhere in JuliaApproximation where this has been done before? Maybe a commit where the necessary changes were made? (perhaps I am too optimistic here). Thank you!

[dlfivefifty]

This is more basic Julia. Eg 0.0 vs zero(T).

some tricks are to use @code_warntype to decide whether the types are inferred (and so code is fast). For unit tests one uses @inferred

[ioannisPApapadopoulos]

Okay! So I think it's working in the sense that giving initial guesses that are BigFloat seems to work. Is there something more subtle that I may have missed?

[ioannisPApapadopoulos]

Changes made! Thanks for going through the code @dlfivefifty

[dlfivefifty]

For the record: we need the following branch

JuliaDiff/ForwardDiff.jl#541

I'll try to get this working with a custom Manifest.toml

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (master@edca52b). Click here to learn what that means.
The diff coverage is n/a.

@@            Coverage Diff            @@
##             master       #5   +/-   ##
=========================================
  Coverage          ?   97.29%           
=========================================
  Files             ?        1           
  Lines             ?       74           
  Branches          ?        0           
=========================================
  Hits              ?       72           
  Misses            ?        2           
  Partials          ?        0           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update edca52b...a71440a. Read the comment docs.