Function space SBP operators

[JoshuaLampert]

This is a very first draft for tackling #132. For the construction of the FSBP an optimization problem needs to be solved. I plan to use Optim.jl for this. Thus, I added a package extension to avoid adding Optim.jl as a direct dependency. In order to discuss some design choices, I already created this PR.

• In the current draft, I use the MatrixDerivativeOperator to represent the functional space operators, i.e. FunctionSpaceOperator just returns MatrixDerivativeOperator. This is easy, but I'm not sure if this is the best design choice on the long run, e.g., this way we "forget" about the basis. A more versatile alternative would be to define a new struct that holds a MatrixDerivativeOperator and potentially more such as the basis for potential later use. In this case, several methods (such as mass_matrix, integrate, *...) could just call the corresponding methods of the MatrixDerivativeOperator.

• Another question concerns the derivatives of the basis functions, which are needed for the construction process. From a user perspective it would be most convenient if they would not need to be specified explicitly. Do you think, we could use ForwardDiff.jl (i.e. the package extension consists of ForwardDiff.jl and Optim.jl) to compute the derivatives or would that mean too much overhead?

I would be interested in your opinion on these two points @ranocha and if the current approach looks reasonable to you.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 90.97%. Comparing base (96e695d) to head (b62f40b).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #268      +/-   ##
==========================================
+ Coverage   90.70%   90.97%   +0.26%     
==========================================
  Files          33       35       +2     
  Lines        5112     5262     +150     
==========================================
+ Hits         4637     4787     +150     
  Misses        475      475              

Flag	Coverage Δ
unittests	90.97% <100.00%> (+0.26%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[coveralls]

coverage: 91.819% (+0.2%) from 91.609%
when pulling 03107f9 on JoshuaLampert:FSBP-extension
into 96e695d on ranocha:main.

[ranocha]

Thanks a lot for the PR!

• I think a MatrixDerivativeOperator is fine for now. Do you have any specific applications using the basis functions later in mind? We could also document that the return value will implement the general interface (like mul! etc.) but leave the concrete type of the return value open - or mention that it's an implementation detail that may change in the future.
• I'm fine using ForwardDiff.jl for this.

[JoshuaLampert]

Thanks for the answers!

• I don't have any specific application in mind, but I thought it might be worth considering the alternatives before one realizes that a more flexible approach would be better. But your suggestion to document that the general interface will be implemented, but the concrete type may change in the future sounds like a good way.
• Ok, I will use ForwardDiff.jl for the derivatives.

[JoshuaLampert]

A first working version is ready. I can (almost) reproduce the results from the paper, e.g. (cf. p. 12),

julia> using SummationByPartsOperators, Optim
Precompiling SummationByPartsOperatorsOptimForwardDiffExt
  1 dependency successfully precompiled in 2 seconds. 120 already precompiled.

julia> D = FunctionSpaceOperator((one, identity, exp), 0.0, 1.0, collect(LinRange(0.0, 1.0, 5)), GlaubitzNordströmÖffner2023())
Matrix-based first-derivative operator {T=Float64} on 5 nodes in [0.0, 1.0]

julia> round.(Matrix(mass_matrix(D)) * D.D, digits = 4)
5×5 Matrix{Float64}:
 -0.5      0.653   -0.035   -0.1927   0.0748
 -0.653    0.0      0.3198   0.5238  -0.1907
  0.035   -0.3198   0.0      0.3215  -0.0367
  0.1927  -0.5238  -0.3215   0.0      0.6526
 -0.0748   0.1907   0.0367  -0.6526   0.5

julia> round.(mass_matrix(D), digits = 4)
5×5 LinearAlgebra.Diagonal{Float64, Vector{Float64}}:
 0.076   ⋅       ⋅       ⋅       ⋅ 
  ⋅     0.3621   ⋅       ⋅       ⋅ 
  ⋅      ⋅      0.1245   ⋅       ⋅ 
  ⋅      ⋅       ⋅      0.3609   ⋅ 
  ⋅      ⋅       ⋅       ⋅      0.0766

Up to now, I have focused on the "make it work" phase by trying to be as close as possible to the mathematical formulation. The "make it fast"(er) phase is yet to come. The optimization process is still taking too long for higher number of nodes, but there is probably a lot of potential for performance improvements. Before I dig deeper into some performance improvement, I would like to see your feedback @ranocha. Does this look ok to you in general or do you think the implementation needs some major alterations?
Up to now I have only added a few tests to cover the most basic. I had a brief look in SBP_operators_test.jl, but was not quite sure, which of these tests also make sense in this case. (I also didn't find any tests like @test M * D.D + D.D' * M ≈ B, which I found a bit surprising as this would have been one of the first things I would like to see tested for the SBP operators 😅)
I have also noticed that the optimization needs a rather small tolerance for the gradient to be able to produce accurate results. I was also wondering about whether we should provide the user some possibility to influence the optimization, e.g., by allowing an additional kwarg options that is passed as Optim.Options to the optimization (with some reasonable default value), or by allowing to use a custom sigmoid function.

[ranocha]

Nice!

Up to now, I have focused on the "make it work" phase by trying to be as close as possible to the mathematical formulation. The "make it fast"(er) phase is yet to come. The optimization process is still taking too long for higher number of nodes, but there is probably a lot of potential for performance improvements. Before I dig deeper into some performance improvement, I would like to see your feedback @ranocha. Does this look ok to you in general or do you think the implementation needs some major alterations?

I think the general approach is fine. I just have some high-level comments right now. The lower-level implementation details could be touched later - we can also discuss this together.

Up to now I have only added a few tests to cover the most basic. I had a brief look in SBP_operators_test.jl, but was not quite sure, which of these tests also make sense in this case. (I also didn't find any tests like @test M * D.D + D.D' * M ≈ B, which I found a bit surprising as this would have been one of the first things I would like to see tested for the SBP operators 😅)

Fair enough - feel free to add tests like this to the other files 😅

I have also noticed that the optimization needs a rather small tolerance for the gradient to be able to produce accurate results. I was also wondering about whether we should provide the user some possibility to influence the optimization, e.g., by allowing an additional kwarg options that is passed as Optim.Options to the optimization (with some reasonable default value), or by allowing to use a custom sigmoid function.

That sounds reasonable.

[JoshuaLampert]

Thanks a lot for your comments!

Fair enough - feel free to add tests like this to the other files 😅

Ok, I can add some tests, but will do that in a future PR.

I have also noticed that the optimization needs a rather small tolerance for the gradient to be able to produce accurate results. I was also wondering about whether we should provide the user some possibility to influence the optimization, e.g., by allowing an additional kwarg options that is passed as Optim.Options to the optimization (with some reasonable default value), or by allowing to use a custom sigmoid function.

That sounds reasonable.

I added options as an additional kwarg (1b8d85e).

[JoshuaLampert]

Since the optimization is doing weird things and does not work in the Julia 1.6 test, we decided to restrict this feature to Julia v1.9, see aa868bd.

[JoshuaLampert]

It is much faster now, but N = 50 nodes still takes ~8.8 seconds on my laptop to create the operators. This scales pretty poorly in N (which is expected since daij_drhok is of size $N \cdot L \cdot K = O(N^3\cdot K)$ ).

[JoshuaLampert]

Do you have any other reasonable tests in mind, I could add, @ranocha?

[JoshuaLampert]

I would also like to stress that the optimization can be very vulnerable. Consider the following example:

julia> using SummationByPartsOperators, Optim

julia> nodes = [0.0, rand(3)..., 1.0]
5-element Vector{Float64}:
 0.0
 0.05651765804974351
 0.29863912793993175
 0.657232875727107
 1.0

julia> basis_functions = [one, identity, exp]
3-element Vector{Function}:
 one (generic function with 45 methods)
 identity (generic function with 1 method)
 exp (generic function with 21 methods)

julia> D = function_space_operator(basis_functions, first(nodes), last(nodes), nodes, GlaubitzNordströmÖffner2023())
Matrix-based first-derivative operator {T=Float64} on 5 nodes in [0.0, 1.0]

julia> D * ones(5)
5-element Vector{Float64}:
 -2038.0570044219494
    -1.073365630488432e-6
     6.732198513992316e-6
     3.664382162149593e-6
    -2.1692864620170837e-5

which is very far from zeros(5). The optimization reports:

 * Status: failure (reached maximum number of iterations)

 * Candidate solution
    Final objective value:     2.677287e-10

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 8.05e-04 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.85e-05 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.90e-22 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 2.58e-12 ≰ 0.0e+00
    |g(x)|                 = 1.17e-11 ≰ 1.0e-14

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    10000
    f(x) calls:    26276
    ∇f(x) calls:   26276

So even with an objective value of ~2.7e-10, which is quite small, the operator is very far from being exact. Increasing the number of iterations significantly and decreasing the tolerance doesn't help:

julia> D = function_space_operator(basis_functions, first(nodes), last(nodes), nodes, GlaubitzNordströmÖffner2023(); options = Optim.Options(g_tol = 1e-16, iterations = 200000))
 * Status: success

 * Candidate solution
    Final objective value:     2.677287e-10

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 2.01e-07 ≰ 0.0e+00
    |x - x'|/|x'|          = 5.40e-09 ≰ 0.0e+00
    |f(x) - f(x')|         = 9.38e-23 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 3.50e-13 ≰ 0.0e+00
    |g(x)|                 = 8.08e-17 ≤ 1.0e-16

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    102609
    f(x) calls:    270140
    ∇f(x) calls:   270140

Matrix-based first-derivative operator {T=Float64} on 5 nodes in [0.0, 1.0]

julia> D * ones(5)
5-element Vector{Float64}:
 -2.46283614165e10
 -1.0733685033015306e-6
  6.7321802603159675e-6
  3.664360717303694e-6
 -2.1692791865035588e-5

I think this is a general issue with the optimization process and thus I think a deeper theoretical insight in the optimization would be useful.
To sum up, if we want to merge this, I would suggest to put a warning in the docstring that this is still an experimental feature and might change in the future.

[ranocha]

To sum up, if we want to merge this, I would suggest to put a warning in the docstring that this is still an experimental feature and might change in the future.

Yes, that's a very good suggestion 👍

[ranocha]

I would also like to stress that the optimization can be very vulnerable.

It looks like this is actually a question for further research and not only an implementation question.

[JoshuaLampert]

It looks like this is actually a question for further research and not only an implementation question.

Yes, definitely. I was just trying to compute the Hessian of the function to check convexity.

[ranocha]

Do you have any other reasonable tests in mind, I could add, @ranocha?

The tests look reasonable 👍

[JoshuaLampert]

To sum up, if we want to merge this, I would suggest to put a warning in the docstring that this is still an experimental feature and might change in the future.

Yes, that's a very good suggestion 👍

Done in 38081b1.

[ranocha]

Thanks!

[JoshuaLampert]

FYI @phioeffn @jglaubitz