Internals

CONSTRUCTION AND/OR APPLICATION OF A SINGLE HOUSEHOLDER TRANSFORMATION Q = I + U*(U**T)/B

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().

CONSTRUCTION AND/OR APPLICATION OF A SINGLE HOUSEHOLDER TRANSFORMATION Q = I + U*(U**T)/B

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

x = nnls(A, b; ...)

Solves non-negative least-squares problem by the active set method of Lawson & Hanson (1974).

Optional arguments:

- max_iter: maximum number of iterations (counts inner loop iterations)

References:

- Lawson, C.L. and R.J. Hanson, Solving Least-Squares Problems
- Prentice-Hall, Chapter 23, p. 161, 1974

COMPUTE ORTHOGONAL ROTATION MATRIX The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 12, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

COMPUTE MATRIX  (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
                (-S,C)         (-S,C)(B)   (   0          )
COMPUTE SIG = SQRT(A**2+B**2)
    SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
    SIG MAY BE IN THE SAME LOCATION AS A OR B .

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 15, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

Algorithm NNLS: NONNEGATIVE LEAST SQUARES

The original version of this code was developed by Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 1973 JUN 15, and published in the book "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. Revised FEB 1995 to accompany reprinting of the book by SIAM.

GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM A * X = B SUBJECT TO X .GE. 0

struct NormalSpline{n, T <: Real, RK <: ReproducingKernel_0} <: AbstractNormalSpline{n,T,RK}

Define a structure containing full information of a normal spline

Fields

• _kernel: a reproducing kernel spline was built with
• _nodes: transformed function value nodes
• _values: function values at interpolation nodes
• _d_nodes: transformed function directional derivative nodes
• _d_dirs: normalized derivative directions
• _d_values: function directional derivative values
• _mu: spline coefficients
• _rhs: right-hand side of the problem gram * mu = rhs
• _gram: Gram matrix of the problem gram * mu = rhs
• _chol: Cholesky factorization of the Gram matrix
• _cond: estimation of the Gram matrix condition number
• _min_bound: minimal bounds of the original node locations area
• _max_bound: maximal bounds of the original node locations area
• _scale: factor of transforming the original node locations into unit hypercube

struct RK_H0{T} <: ReproducingKernel_0

Defines a type of reproducing kernel of Bessel Potential space $H^{n/2 + 1/2}_ε (R^n)$ ('Basic Matérn kernel'):

\[V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|) \, .\]

Fields

• ε::T: 'scaling parameter' from the Bessel Potential space definition, it may be omitted in the struct constructor otherwise it must be greater than zero

struct RK_H1{T} <: ReproducingKernel_1

Defines a type of reproducing kernel of Bessel Potential space $H^{n/2 + 3/2}_ε (R^n)$ ('Linear Matérn kernel'):

\[V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|) (1 + \varepsilon |\xi - \eta|) \, .\]

Fields

• ε::T: 'scaling parameter' from the Bessel Potential space definition, it may be omitted in the struct constructor otherwise it must be greater than zero

struct RK_H2{T} <: ReproducingKernel_2

Defines a type of reproducing kernel of Bessel Potential space $H^{n/2 + 5/2}_ε (R^n)$ ('Quadratic Matérn kernel'):

\[V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|) (3 + 3\varepsilon |\xi - \eta| + \varepsilon ^2 |\xi - \eta| ^2 ) \, .\]

Fields

• ε::T: 'scaling parameter' from the Bessel Potential space definition, it may be omitted in the struct constructor otherwise it must be greater than zero

Get estimation of the Gram matrix condition number Brás, C.P., Hager, W.W. & Júdice, J.J. An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP 20, 791–809 (2012). https://link.springer.com/article/10.1007/s11750-010-0161-9

construct(spline::AbstractNormalSpline{n,T,RK}, values::AbstractVector{T}, d_values::AbstractVector{T}) where {n, T <: Real, RK <: ReproducingKernel_1}

Construct the spline by calculating its coefficients and completely initializing the NormalSpline object.

Arguments

• spline: the partly initialized NormalSpline object returned by prepare function.
• values: function values at nodes nodes.
• d_values: function directional derivative values at d_nodes nodes.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function to interpolate the data to required points.

construct(spline::AbstractNormalSpline{n,T,RK}, values::AbstractVector{T}) where {n, T <: Real, RK <: ReproducingKernel_0}

Construct the spline by calculating its coefficients and completely initializing the NormalSpline object.

Arguments

• spline: the partly initialized NormalSpline object returned by prepare function.
• values: function values at nodes nodes.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function to interpolate the data to required points.

estimate_accuracy(spline::AbstractNormalSpline{n,T,RK}) where {n, T <: Real, RK <: ReproducingKernel_0}

Assess accuracy of interpolation results by analyzing residuals.

Arguments

• spline: the NormalSpline object returned by construct or interpolate function.

Returns

• An estimation of the number of significant digits in the interpolation result.

estimate_cond(spline::AbstractNormalSpline{n,T,RK}) where {n, T <: Real, RK <: ReproducingKernel_0}

Get an estimation of the Gram matrix condition number. It needs the spline object is prepared and requires O(N^2) operations. (C. Brás, W. Hager, J. Júdice, An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP Vol.20, No.3, 2012.)

Arguments

• spline: the NormalSpline object returned by prepare, construct or interpolate function.

Returns

• An estimation of the Gram matrix condition number.

estimate_epsilon(nodes::AbstractMatrix{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Get the estimation of the 'scaling parameter' of Bessel Potential space the spline being built in. It coincides with the result returned by get_epsilon function.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• kernel: reproducing kernel of Bessel potential space the normal spline will be constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• Estimation of ε.

estimate_epsilon(nodes::AbstractVector{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Get an the estimation of the 'scaling parameter' of Bessel Potential space the 1D spline is being built in. It coincides with the result returned by get_epsilon function.

Arguments

• nodes: The function value nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• Estimation of ε.

estimate_epsilon(nodes::AbstractMatrix{T}, d_nodes::AbstractMatrix{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Get an the estimation of the 'scaling parameter' of Bessel Potential space the spline being built in. It coincides with the result returned by get_epsilon function.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• d_nodes: The function directional derivative nodes. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline will be constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• Estimation of ε.

estimate_epsilon(nodes::AbstractVector{T}, d_nodes::AbstractVector{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Get an the estimation of the 'scaling parameter' of Bessel Potential space the 1D spline is being built in. It coincides with the result returned by get_epsilon function.

Arguments

• nodes: The function value nodes.
• d_nodes: The function derivative nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• Estimation of ε.

evaluate(spline::AbstractNormalSpline{1,T,RK}, point::T) where {T <: Real, RK <: ReproducingKernel_0}

Evaluate the 1D spline value at the point location.

Arguments

• spline: the NormalSpline object returned by interpolate or construct function.
• point: location at which spline value is evaluating.

Returns

• Spline value at the point location.

evaluate(spline::AbstractNormalSpline{n,T,RK}, points::AbstractMatrix{T}) where {n, T <: Real, RK <: ReproducingKernel_0}

Evaluate the spline values at the locations defined in points.

Arguments

• spline: theNormalSplineobject returned byinterpolateorconstruct` function.
• points: locations at which spline values are evaluating. This should be an n×m matrix, where n is dimension of the sampled space and m is the number of locations where spline values are evaluating. It means that each column in the matrix defines one location.

Returns

• Vector{T} of the spline values at the locations defined in points.

evaluate(spline::AbstractNormalSpline{n,T,RK}, points::AbstractVector{T}) where {n, T <: Real, RK <: ReproducingKernel_0}

Evaluate the 1D spline values/value at the points locations.

Arguments

• spline: the NormalSpline object returned by interpolate or construct function.
• points: locations at which spline values are evaluating. This should be a vector of size m where m is the number of evaluating points.

Returns

• Spline value at the point location.

evaluate_derivative(spline::AbstractNormalSpline{1,T,RK}, point::T) where {T <: Real, RK <: ReproducingKernel_0}

Evaluate the 1D spline derivative at the point location.

Arguments

• spline: the NormalSpline object returned by interpolate or construct function.
• point: location at which spline derivative is evaluating.

Note: Derivative of spline built with reproducing kernel RK_H0 does not exist at the spline nodes.

Returns

• The spline derivative value at the point location.

evaluate_gradient(spline::AbstractNormalSpline{n,T,RK}, point::AbstractVector{T}) where {n, T <: Real, RK <: ReproducingKernel_0}

Evaluate gradient of the spline at the location defined in point.

Arguments

• spline: the NormalSpline object returned by interpolate or construct function.
• point: location at which gradient value is evaluating. This should be a vector of size n, where n is dimension of the sampled space.

Note: Gradient of spline built with reproducing kernel RK_H0 does not exist at the spline nodes.

Returns

• Vector{T} - gradient of the spline at the location defined in point.

evaluate_one(spline::AbstractNormalSpline{n,T,RK}, point::AbstractVector{T}) where {n, T <: Real, RK <: ReproducingKernel_0}

Evaluate the spline value at the point location.

Arguments

• spline: the NormalSpline object returned by interpolate or construct function.
• point: location at which spline value is evaluating. This should be a vector of size n, where n is dimension of the sampled space.

Returns

• The spline value at the location defined in point.

get_cond(nodes::AbstractMatrix{T}, d_nodes::AbstractMatrix{T}, d_dirs::AbstractMatrix{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Get a value of the Gram matrix spectral condition number. It is obtained by means of the matrix SVD decomposition and requires $O(N^3)$ operations.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• d_nodes: The function directional derivatives nodes. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes.
• d_dirs: Directions of the function directional derivatives. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes. It means that each column in the matrix defines one direction of the function directional derivative.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• A value of the Gram matrix spectral condition number.

get_cond(nodes::AbstractMatrix{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Get a value of the Gram matrix spectral condition number. It is obtained by means of the matrix SVD decomposition and requires $O(N^3)$ operations.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• A value of the Gram matrix spectral condition number.

get_epsilon(spline::AbstractNormalSpline{n,T,RK}) where {n, T <: Real, RK <: ReproducingKernel_0}

Get the 'scaling parameter' of Bessel Potential space the spline was built in.

Arguments

• spline: the NormalSpline object returned by prepare, construct or interpolate function.

Returns

• The 'scaling parameter' ε.

interpolate(nodes::AbstractVector{T}, values::AbstractVector{T}, d_nodes::AbstractVector{T}, d_values::AbstractVector{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Prepare and construct the 1D interpolating normal spline.

Arguments

• nodes: function value interpolation nodes. This should be an n₁ vector where n₁ is the number of function value nodes.
• values: function values at nodes nodes.
• d_nodes: The function derivatives nodes. This should be an n₂ vector where n₂ is the number of function derivatives nodes.
• d_values: function derivative values at d_nodes nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function.

interpolate(nodes::AbstractMatrix{T}, values::AbstractVector{T}, d_nodes::AbstractMatrix{T}, d_dirs::AbstractMatrix{T}, d_values::AbstractVector{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Prepare and construct the spline.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• values: function values at nodes nodes.
• d_nodes: The function directional derivative nodes. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes.
• d_dirs: Directions of the function directional derivatives. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes. It means that each column in the matrix defines one direction of the function directional derivative.
• d_values: function directional derivative values at d_nodes nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function.

interpolate(nodes::AbstractMatrix{T}, values::AbstractVector{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Prepare and construct the spline.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• values: function values at nodes nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function.

interpolate(nodes::AbstractVector{T}, values::AbstractVector{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Prepare and construct the 1D spline.

Arguments

• nodes: function value interpolation nodes. This should be an n₁ vector where n₁ is the number of function value nodes.
• values: function values at n₁ interpolation nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The completely initialized NormalSpline object that can be passed to evaluate function.

prepare(nodes::AbstractMatrix{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Prepare the spline by constructing and factoring a Gram matrix of the interpolation problem. Initialize the NormalSpline object.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The partly initialized NormalSpline object that must be passed to construct function in order to complete the spline initialization.

prepare(nodes::AbstractVector{T}, kernel::RK = RK_H0()) where {T <: Real, RK <: ReproducingKernel_0}

Prepare the 1D spline by constructing and factoring a Gram matrix of the interpolation problem. Initialize the NormalSpline object.

Arguments

• nodes: function value interpolation nodes. This should be an n₁ vector where n₁ is the number of function value nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H0 if the spline is constructing as a continuous function, RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The partly initialized NormalSpline object that must be passed to construct function in order to complete the spline initialization.

prepare(nodes::AbstractMatrix{T}, d_nodes::AbstractMatrix{T}, d_dirs::AbstractMatrix{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Prepare the spline by constructing and factoring a Gram matrix of the interpolation problem. Initialize the NormalSpline object.

Arguments

• nodes: The function value nodes. This should be an n×n_1 matrix, where n is dimension of the sampled space and n₁ is the number of function value nodes. It means that each column in the matrix defines one node.
• d_nodes: The function directional derivatives nodes. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes.
• d_dirs: Directions of the function directional derivatives. This should be an n×n_2 matrix, where n is dimension of the sampled space and n₂ is the number of function directional derivative nodes. It means that each column in the matrix defines one direction of the function directional derivative.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The partly initialized NormalSpline object that must be passed to construct function in order to complete the spline initialization.

prepare(nodes::AbstractVector{T}, d_nodes::AbstractVector{T}, kernel::RK = RK_H1()) where {T <: Real, RK <: ReproducingKernel_1}

Prepare the 1D interpolating normal spline by constructing and factoring a Gram matrix of the problem. Initialize the NormalSpline object.

Arguments

• nodes: function value interpolation nodes. This should be an n₁ vector where n₁ is the number of function value nodes.
• d_nodes: The function derivatives nodes. This should be an n₂ vector where n₂ is the number of function derivatives nodes.
• kernel: reproducing kernel of Bessel potential space the normal spline is constructed in. It must be a struct object of the following type: RK_H1 if the spline is constructing as a differentiable function, RK_H2 if the spline is constructing as a twice differentiable function.

Returns

• The partly initialized NormalSpline object that must be passed to construct function in order to complete the spline initialization.

LinearAlgebra.cholesky!(C::ElasticCholesky, v::AbstractVector{T}) where {T}

Update the Cholesky factorization C as if the column v (and by symmetry, the corresponding row vᵀ) were inserted into the underlying matrix A. Specifically, let L be the lower-triangular cholesky factor of A such that A = LLᵀ, and let v = [d; γ] such that the new matrix A⁺ is given by

A⁺ = [A  d]
     [dᵀ γ].

Then, the corresponding updated cholesky factor L⁺ of ⁺ is:

L⁺ = [L  0]
     [eᵀ α]

where e = L⁻¹d, α = √τ, and τ = γ - e⋅e > 0. If τ ≤ 0, then A⁺ is not positive definite.

See: https://igorkohan.github.io/NormalHermiteSplines.jl/dev/Normal-Splines-Method/#Algorithms-for-updating-Cholesky-factorization