Add least squares approximation (#93)

* allow using different node set for centers

* add example and test

* add unit tests

examples/interpolation/least_squares_2d.jl

@@ -0,0 +1,28 @@
+using KernelInterpolation
+using Plots
+
+# interpolate Franke function
+function f(x)
+    0.75 * exp(-0.25 * ((9 * x[1] - 2)^2 + (9 * x[2] - 2)^2)) +
+    0.75 * exp(-(9 * x[1] + 1)^2 / 49 - (9 * x[2] + 1) / 10) +
+    0.5 * exp(-0.25 * ((9 * x[1] - 7)^2 + (9 * x[2] - 3)^2)) -
+    0.2 * exp(-(9 * x[1] - 4)^2 - (9 * x[2] - 7)^2)
+end
+
+n = 1089
+nodeset = random_hypercube(n; dim = 2)
+values = f.(nodeset) .+ 0.03 * randn(n)
+M = 81
+centers = random_hypercube(M; dim = 2)
+
+kernel = ThinPlateSplineKernel{dim(nodeset)}()
+ls = interpolate(nodeset, centers, values, kernel)
+itp = interpolate(nodeset, values, kernel)
+
+N = 40
+many_nodes = homogeneous_hypercube(N; dim = 2)
+
+p1 = plot(many_nodes, ls, st = :surface, training_nodes = false)
+p2 = plot(many_nodes, itp, st = :surface, training_nodes = false)
+
+plot(p1, p2, layout = (2, 1))

src/discretization.jl

@@ -53,7 +53,8 @@ function solve_stationary(spatial_discretization::SpatialDiscretization{Dim, Rea
     # Do not support additional polynomial basis for now
     xx = polyvars(Dim)
     ps = monomials(xx, 0:-1)
-    return Interpolation(kernel, merge(nodeset_inner, nodeset_boundary), c, system_matrix,
+    nodeset = merge(nodeset_inner, nodeset_boundary)
+    return Interpolation(kernel, nodeset, nodeset, c, system_matrix,
                          ps, xx)
 end
 

src/interpolation.jl

@@ -19,6 +19,7 @@ struct Interpolation{Kernel, Dim, RealT, A, Monomials, PolyVars} <:
        AbstractInterpolation{Kernel, Dim, RealT}
     kernel::Kernel
     nodeset::NodeSet{Dim, RealT}
+    centers::NodeSet{Dim, RealT}
     c::Vector{RealT}
     system_matrix::A
     ps::Monomials
@@ -69,7 +70,7 @@ interpolant.
 
 See also [`coefficients`](@ref) and [`polynomial_coefficients`](@ref).
 """
-kernel_coefficients(itp::Interpolation) = itp.c[eachindex(nodeset(itp))]
+kernel_coefficients(itp::Interpolation) = itp.c[eachindex(itp.centers)]
 
 """
     polynomial_coefficients(itp::Interpolation)
@@ -79,7 +80,7 @@ interpolant.
 
 See also [`coefficients`](@ref) and [`kernel_coefficients`](@ref).
 """
-polynomial_coefficients(itp::Interpolation) = itp.c[(length(nodeset(itp)) + 1):end]
+polynomial_coefficients(itp::Interpolation) = itp.c[(length(itp.centers) + 1):end]
 
 """
     polynomial_basis(itp::Interpolation)
@@ -121,7 +122,7 @@ or collocation is used.
 system_matrix(itp::Interpolation) = itp.system_matrix
 
 @doc raw"""
-    interpolate(nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel))
+    interpolate(nodeset, centers = nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel))
 
 Interpolate the `values` evaluated at the nodes in the `nodeset` to a function using the kernel `kernel`
 and polynomials up to a order `m` (i.e. degree - 1), i.e., determine the coefficients ``c_j`` and ``d_k`` in the expansion
@@ -135,8 +136,10 @@ maximum degree of `m - 1`. If `m = 0`, no polynomial is added. The additional co
     \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q = \begin{pmatrix}m - 1 + d\\d\end{pmatrix}
 ```
 are enforced. Returns an [`Interpolation`](@ref) object.
+If `centers` is provided, the interpolant is a least squares approximation with the centers used for the basis.
 """
-function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT},
+function interpolate(nodeset::NodeSet{Dim, RealT},
+                     values::Vector{RealT},
                      kernel = GaussKernel{Dim}();
                      m = order(kernel)) where {Dim, RealT}
     @assert dim(kernel) == Dim
@@ -149,24 +152,46 @@ function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT},
     k_matrix = kernel_matrix(nodeset, kernel)
     p_matrix = polynomial_matrix(nodeset, ps)
     system_matrix = [k_matrix p_matrix
-                     transpose(p_matrix) zeros(q, q)]
+                     p_matrix' zeros(q, q)]
     b = [values; zeros(q)]
     symmetric_system_matrix = Symmetric(system_matrix)
     c = symmetric_system_matrix \ b
-    return Interpolation(kernel, nodeset, c, symmetric_system_matrix, ps, xx)
+    return Interpolation(kernel, nodeset, nodeset, c, symmetric_system_matrix, ps, xx)
+end
+
+# Least squares approximation
+function interpolate(nodeset::NodeSet{Dim, RealT}, centers::NodeSet{Dim, RealT},
+                     values::Vector{RealT},
+                     kernel = GaussKernel{Dim}();
+                     m = order(kernel)) where {Dim, RealT}
+    @assert dim(kernel) == Dim
+    n = length(nodeset)
+    @assert length(values) == n
+    xx = polyvars(Dim)
+    ps = monomials(xx, 0:(m - 1))
+    q = length(ps)
+
+    k_matrix = kernel_matrix(nodeset, centers, kernel)
+    p_matrix1 = polynomial_matrix(nodeset, ps)
+    p_matrix2 = polynomial_matrix(centers, ps)
+    system_matrix = [k_matrix p_matrix1
+                     p_matrix2' zeros(q, q)]
+    b = [values; zeros(q)]
+    c = system_matrix \ b
+    return Interpolation(kernel, nodeset, centers, c, system_matrix, ps, xx)
 end
 
 # Evaluate interpolant
 function (itp::Interpolation)(x)
     s = 0
     kernel = interpolation_kernel(itp)
-    xs = nodeset(itp)
+    xis = itp.centers
     c = kernel_coefficients(itp)
     d = polynomial_coefficients(itp)
     ps = polynomial_basis(itp)
     xx = polyvars(itp)
     for j in eachindex(c)
-        s += c[j] * kernel(x, xs[j])
+        s += c[j] * kernel(x, xis[j])
     end
 
     for k in eachindex(d)
@@ -184,22 +209,22 @@ end
 function (diff_op_or_pde::Union{AbstractDifferentialOperator, AbstractStationaryEquation})(itp::Interpolation,
                                                                                            x)
     kernel = interpolation_kernel(itp)
-    xs = nodeset(itp)
+    xis = itp.centers
     c = kernel_coefficients(itp)
     s = zero(eltype(x))
     for j in eachindex(c)
-        s += c[j] * diff_op_or_pde(kernel, x, xs[j])
+        s += c[j] * diff_op_or_pde(kernel, x, xis[j])
     end
     return s
 end
 
 function (g::Gradient)(itp::Interpolation, x)
     kernel = interpolation_kernel(itp)
-    xs = nodeset(itp)
+    xis = itp.centers
     c = kernel_coefficients(itp)
     s = zero(x)
     for j in eachindex(c)
-        s += c[j] * g(kernel, x, xs[j])
+        s += c[j] * g(kernel, x, xis[j])
     end
     return s
 end
@@ -223,8 +248,8 @@ function kernel_inner_product(itp1, itp2)
     @assert kernel == interpolation_kernel(itp2)
     c_f = kernel_coefficients(itp1)
     c_g = kernel_coefficients(itp2)
-    xs = nodeset(itp1)
-    xis = nodeset(itp2)
+    xs = itp1.centers
+    xis = itp2.centers
     s = 0
     for i in eachindex(c_f)
         for j in eachindex(c_g)
@@ -274,7 +299,8 @@ function (titp::TemporalInterpolation)(t)
     # Do not support additional polynomial basis for now
     xx = polyvars(dim(semi))
     ps = monomials(xx, 0:-1)
-    itp = Interpolation(kernel, merge(nodeset_inner, nodeset_boundary), c,
+    nodeset = merge(nodeset_inner, nodeset_boundary)
+    itp = Interpolation(kernel, nodeset, nodeset, c,
                         semi.cache.mass_matrix, ps, xx)
     return itp
 end

src/visualization.jl

@@ -106,13 +106,13 @@ end
     elseif dim(nodeset) == 2
         if training_nodes
             @series begin
-                x = values_along_dim(itp.nodeset, 1)
-                y = values_along_dim(itp.nodeset, 2)
+                x = values_along_dim(itp.centers, 1)
+                y = values_along_dim(itp.centers, 2)
                 seriestype := :scatter
                 markershape --> :star
                 markersize --> 10
                 label --> "training nodes"
-                x, y, itp.(itp.nodeset)
+                x, y, itp.(itp.centers)
             end
         end
         @series begin

test/test_examples_interpolation.jl

@@ -118,3 +118,14 @@ end
     @test_include_example(joinpath(EXAMPLES_DIR, "interpolation_2d_condition.jl"),
                           ns=5:10)
 end
+
+@testitem "least_squares_2d.jl" setup=[
+    Setup,
+    AdditionalImports,
+    InterpolationExamples
+] begin
+    @test_include_example(joinpath(EXAMPLES_DIR, "least_squares_2d.jl"),
+                          l2=1.2759520194191292, linf=0.19486087346749836,
+                          l2_ls=0.5375130503454387, linf_ls=0.06810374254243684,
+                          interpolation_test=false, least_square_test=true)
+end

test/test_unit.jl

@@ -675,9 +675,37 @@ end
     @test length(polynomial_basis(itp)) ==
           binomial(order(itp) - 1 + dim(nodes), dim(nodes))
     @test system_matrix(itp) isa Symmetric
+    @test size(system_matrix(itp)) == (7, 7)
     @test isapprox(itp([0.5, 0.5]), 1.0)
     @test isapprox(kernel_norm(itp), 0.0)
 
+    # Least squares approximation
+    centers = NodeSet([0.0 0.0
+                       1.0 0.0
+                       0.0 1.0])
+    itp = @test_nowarn interpolate(nodes, centers, ff, kernel)
+    expected_coefficients = [
+        0.0,
+        0.0,
+        0.0,
+        0.0,
+        1.0,
+        1.0]
+    coeffs = coefficients(itp)
+    @test length(coeffs) == length(expected_coefficients)
+    for i in eachindex(coeffs)
+        @test isapprox(coeffs[i], expected_coefficients[i], atol = 1e-15)
+    end
+    @test order(itp) == order(kernel)
+    @test length(kernel_coefficients(itp)) == length(itp.centers)
+    @test length(polynomial_coefficients(itp)) == order(itp) + 1
+    @test length(polynomial_basis(itp)) ==
+          binomial(order(itp) - 1 + dim(nodes), dim(nodes))
+    @test system_matrix(itp) isa Matrix
+    @test size(system_matrix(itp)) == (7, 6)
+    @test isapprox(itp([0.5, 0.5]), 1.0)
+    @test isapprox(kernel_norm(itp), 0.0, atol = 1e-15)
+
     # 1D interpolation and evaluation
     nodes = NodeSet(LinRange(0.0, 1.0, 10))
     f(x) = sinpi(x[1])

test/test_util.jl

@@ -1,21 +1,27 @@
 """
     test_include_example(example; l2=nothing, linf=nothing,
+                         l2_ls=nothing, linf_ls=nothing,
                          atol=1e-12, rtol=sqrt(eps()),
-                         args...)
+                         kwargs...)
 
-Test by calling `include_example(example; parameters...)`.
+Test by calling `include_example(example; kwargs...)`.
 By default, only the absence of error output is checked.
 """
 macro test_include_example(example, args...)
     local l2 = get_kwarg(args, :l2, nothing)
     local linf = get_kwarg(args, :linf, nothing)
+    local l2_ls = get_kwarg(args, :l2_ls, nothing)
+    local linf_ls = get_kwarg(args, :linf_ls, nothing)
     local atol = get_kwarg(args, :atol, 1e-12)
     local rtol = get_kwarg(args, :rtol, sqrt(eps()))
+    local interpolation_test = get_kwarg(args, :interpolation_test, true)
+    local least_square_test = get_kwarg(args, :least_square_test, false)
     local pde_test = get_kwarg(args, :pde_test, false)
     local kwargs = Pair{Symbol, Any}[]
     for arg in args
         if (arg.head == :(=) &&
-            !(arg.args[1] in (:l2, :linf, :atol, :rtol, :pde_test)))
+            !(arg.args[1] in (:l2, :linf, :l2_ls, :linf_ls, :atol, :rtol,
+                              :interpolation_test, :least_square_test, :pde_test)))
             push!(kwargs, Pair(arg.args...))
         end
     end
@@ -32,14 +38,23 @@ macro test_include_example(example, args...)
                 # assumes `many_nodes` and `values` are defined in the example
                 values_test = itp.(nodeset)
                 # Check interpolation at interpolation nodes
-                @test isapprox(norm(values .- values_test, Inf), 0;
-                               atol = $atol, rtol = $rtol)
+                if $interpolation_test
+                    @test isapprox(norm(values .- values_test, Inf), 0;
+                                   atol = $atol, rtol = $rtol)
+                end
                 many_values = f.(many_nodes)
                 many_values_test = itp.(many_nodes)
                 @test isapprox(norm(many_values .- many_values_test), $l2;
                                atol = $atol, rtol = $rtol)
                 @test isapprox(norm(many_values .- many_values_test, Inf), $linf;
                                atol = $atol, rtol = $rtol)
+                if $least_square_test
+                    many_values_ls = ls.(many_nodes)
+                    @test isapprox(norm(many_values .- many_values_ls), $l2_ls;
+                                   atol = $atol, rtol = $rtol)
+                    @test isapprox(norm(many_values .- many_values_ls, Inf), $linf_ls;
+                                   atol = $atol, rtol = $rtol)
+                end
             else
                 # PDE test
                 # assumes `many_nodes`, `nodes_inner` and `nodeset_boundary` are defined in the example