Guide

Installation

The package can be installed with the Julia package manager. From the Julia REPL, type ] to enter the Pkg REPL mode and run

pkg> add VP4Optim

Model definition

VARPRO can be applied to any data model that can be written in the form $\bm{y}\approx \bm{A}\,\bm{c}$. Apparently, the nontrivial part of the model definition depends on how the matrix $\bm{A}\left(\bm{p}\right)$ is actually defined as a function of tunable parameters $\bm{p}$. With respect to optimization though, we sometimes may only want some subset $\bm{x}\subseteq\bm{p}$ to be variable while keeping the rest fixed^[1]. In addition, we want to remain flexible with respect to selecting different subsets of $\bm{p}$ as variable.

To address these requirements, any specific VARPRO model in VP4Optim shall be defined as some subtype of the abstract type VP4Optim.Model{Ny,Nx,Nc,T}

VP4Optim.Model — Type

Model{Ny,Nx,Nc,T}

Abstract supertype of any model specification.

Type parameters

Ny::Int: length of data vector y.
Nx::Int: number of variable parameters.
Nc::Int: number of linear coefficients.
T <: Union{Float64, ComplexF64} is equal to eltype(y)

Some fields are mandatory in any model definition, as shown in the following example of a complex valued model^[2]:

import VP4Optim as VP
using StaticArrays

mutable struct SpecificModel{Ny,Nx,Nc} <: VP.Model{Ny,Nx,Nc,ComplexF64}
    # mandatory fields of any Model instance

    # names of all parameters, variable and fixed 
    # Example: sym == [:a, :b, :c, :d, :e]
    
    sym::Vector{Symbol}     

    # subset of variable parameter names
    # Example: x_sym == [:d, :a, :c] (order defined by x_ind below)
    # Note that this also implies Nx == 3
    
    x_sym::Vector{Symbol}   
    
    # complementary subset of fixed parameter names
    # Example: par_sym == [:e, :b] (order defined by par_ind below)
    
    par_sym::Vector{Symbol} 
    
    # values of all parameters in the order specified by sym
    # Example: val == [a, b, c, d, e]
    
    val::Vector{Float64}
    
    # indices of variable parameters in field val
    # Example: x_ind == SVector{Nx}([4, 1, 3]) (according to the definition of x_sym above)
    # such that
    # x_sym[1] == sym[x_ind[1]] == sym[4] == :d
    # val[x_ind[1]] == val[4] == d
    # x_sym[2] == sym[x_ind[2]] == sym[1] == :a
    # val[x_ind[2]] == val[1] == a
    # x_sym[3] == sym[x_ind[3]] == sym[3] == :c
    # val[x_ind[3]] == val[3] == c
    
    x_ind::SVector{Nx,Int}

    # indices of fixed parameters in field val
    # Example: par_ind == [5, 2] (according to the definition of par_sym above)
    # such that
    # par_sym[1] == sym[par_ind[1]] == sym[5] == :e
    # val[par_ind[1]] == val[5] == e
    # par_sym[2] == sym[par_ind[2]] == sym[2] == :b
    # val[par_ind[2]] == val[2] == b
    
    par_ind::Vector{Int}
    
    # actual data vector
    
    y::SVector{Ny,ComplexF64}            
    
    # == real(y' * y) (by default automatically calculated in generic method VP.y!())
    
    y2::Float64                 

    # optional model-specific information, needed for the constructor
    X::Float64
    Y::Symbol
    time_points::Vector{Float64}
    Z::Vector{ComplexF64}

    # and/or allocated workspace
    # ...
end

Note

For general models, one can include the matrix A::SMatrix{Ny, Nc} as an additional field in any model implementation, since this prevents unnecessary recomputations of A. With this approach, one would let the methods x_changed! and par_changed! trigger updates of the field A.

For proper functioning of the package, neither field A or method A are mandatory though, as long as the methods Bb!, ∂Bb! and ∂∂Bb! are implemented properly (cf. the method descriptions for more details). Depending on the model, this route can often be preferrable to improve numerical performance.

Constructor parameters

Depending on the model, the number of arguments for the Constructor can be large. For this reason, we collect these arguments as fields in a subtype of

VP4Optim.ModPar — Type

ModPar{T}

Abstract supertype of model parameters.

Remarks

Intended to be implemented as a non-mutable structure.

For the example model above, the definition could look like

struct SpeModPar <: VP.ModPar{SpecificModel}
    sym::Vector{Symbol}
    x_sym::Vector{Symbol}
    X::Float64
    Y::Symbol
    time_points::Vector{Float64}
end

Note

Instances of ModPar are only used as arguments for the Constructor.

Internal parameters, like a,b,c,d,e in the example above, are typically set via par! or x! and therefore usually not included in subtypes of ModPar.

For modpar to work, any subtype of ModPar should provide a default constructor without arguments

VP4Optim.ModPar — Method

ModPar(::Type{<: Model})

Construct default parameters for model of type T.

Remark

Must be implemented for each specific model.

like this

function VP.ModPar(::Type{SpecificModel})
    sym_ = [:a, :b, :c, :d, :e]
    x_sym_ = deepcopy(sym)
    X_ = 0.0
    Y_ = :42
    time_points_ = Float64[]
    
    SpeModPar(sym_, x_sym_, X_, Y_, time_points_)
end

Note

There is no need to define ModPar mutable, since the modpar routines are supplied.

To specify model settings, the modpar routines can be used

VP4Optim.modpar — Function

modpar(T::Type{<: Model}; kwargs...)

Arguments

T::Type{<: Model}: Type of the model to be constructed.
kwargs: Arbitrary number of keyword arguments.

Returns

Instance of type <: ModPar{T}.

Remarks

First generates default parameters with ModPar (which must be implemented for each model)

modpar(mp::ModPar{T}; kwargs...) where {T <: Model}

Arguments

mp::ModPar{T}: Existing model parameters to be changed.
kwargs: Arbitrary number of keyword arguments.

Returns

Instance of type <: ModPar{T}.

Remarks

Note that the argument mp is not modified.

which are applied like this

# Generate an instance of ModPar either with default settings ...
smp = modpar(SpecificModel)  # equivalent to smp = VP.ModPar(SpecificModel)
# ... or specific settings via one or more keyword arguments
smp = modpar(SpecificModel; x_sym = [:a, :d], Y_ = :43)

# Parameters of an instance smp can also be changed
smp = modpar(smp; x_sym = [:a, :b, :c], X_ = 1.0, time_points_ = [0, 1, 2])

Note

modpar does not change the supplied argument smp but returns a new one, which must be caught.

Constructor

For given model parameters, a model instance can always be generated with the function

VP4Optim.create_model — Function

create_model(mp::ModPar{T}) where {T <: Model}

Create a new instance of model T with model parameters mp

Remarks

Calls check to guarantee consistency of parameters.
Then calls the constructor T(mp), which must be implemented for each model.

like this

# Generate an instance of SpecificModel
mod = VP.create_model(smp)

Note

For this to work, each model SpecificModel must provide a constuctor SpecificModel(::SpeModPar).
To benefit from the increased performance of StaticArrays.jl, the parameters Ny, Nx, Nc need to be known at compile time. This can be accomplished with the following exemplary setup:

function SpecificModel(smp::SpeModPar)
    # Number of variable parameters
    Nx = length(smp.x_sym)

    # In our example, the number of data points must be equal to the number of time points:
    Ny = length(smp.time_points)

    # If not already fixed by SpecificModel, Nc must also be inferred from smp:
    Nc = some_function_of(smp)

    # Now we can call a constructor, where Ny, Nx, Nc are converted into types
    SpecificModel(Val(Ny), Val(Nx), Val(Nc), smp)
end

function SpecificModel(::Val{Ny}, ::Val{Nx}, ::Val{Nc}, smp) where {Ny, Nx, Nc}
    # name the parameters as desired
    sym = deepcopy(smp.sym)
    
    # specify the variable subset of the parameters
    x_sym = deepcopy(smp.x_sym)
    
    # The remaining mandatory fields can be calculated as follows
    val = zeros(length(sym))
    x_ind = SVector{Nx,Int}(findfirst(x -> x == x_sym[i], sym) for i in 1:Nx)
    par_ind = filter(x -> x ∉ x_ind, 1:length(sym))
    par_sym = sym[par_ind]
    y = SVector{Ny,ComplexF64}(zeros(ComplexF64, Ny))
    y2 = 0.0
    
    # Optionally, initialize further fields, which appear in SpecificModel
    X = smp.X
    Y = smp.Y
    ts = deepcopy(smp.time_points)
    Z = exp.(im * smp.X * ts)

    # Finally, generate an instance of SpecificModel with the default constructor
    SpecificModel{Ny, Nx, Nc}(sym, x_sym, par_sym, val, x_ind, par_ind, y, y2, X, Y, Z, ts)
end

Note

Often, the number of linear coefficients Nc is already determined by the SpecificModel type. In such cases, Nc will not appear as a type parameter.

An example is the biexponential (Nc == 2) decay model in BiExpDecay.jl

mutable struct BiExpDecay{Ny, Nx} <: VP.Model{Ny, Nx, 2, ComplexF64}
...
end

Also the constructors will then need to be adapted accordingly, as shown there.

Note

The method create_model calls the function check before actually calling the constructor.
If model parameters have to fulfil certain requirements, these model-dependent constistency checks should therefore be placed in check:

VP4Optim.check — Function

check(::ModPar)

Throws an Exception, if the supplied constructor parameters are inconsistent.

Arguments

::ModPar: Instance of model parameters to be checked.

Default

Does nothing.

Remarks

Is automatically called by create_model.

In our case, an implementation of check could looks like this

function check(smp::SpeModPar)
    @assert length(smp.sym) == 5
    @assert all(sy -> sy ∈ smp.sym, smp.x_sym)
    @assert smp.X ≥ 0
    # ...
end

Model parameters

Based upon the mandatory fields in any Model instance, the following routines allow access and (partly) modify the fields (after initialization in the constructor).

VP4Optim.N_data — Function

N_data(::Model{Ny,Nx,Nc,T}) where {Ny,Nx,Nc,T}

Returns Ny (number of data)

VP4Optim.N_var — Function

N_var(::Model{Ny,Nx,Nc,T}) where {Ny,Nx,Nc,T}

Returns Nx (number of variable parameters)

VP4Optim.N_coeff — Function

N_coeff(::Model{Ny,Nx,Nc,T}) where {Ny,Nx,Nc,T}

Returns Nc (number of linear coefficients)

VP4Optim.data_type — Function

data_type(::Model{Ny,Nx,Nc,T}) where {Ny,Nx,Nc,T}

Returns T (data type of model)

VP4Optim.sym — Function

sym(mod::Model)

Return all model parameter names (type: Symbol), variable and fixed.

Default

Returns mod.sym::Vector{Symbol}.

VP4Optim.val — Function

val(mod::Model)

Returns all model parameters.

Default

Returns mod.val::Vector{Float64}.

VP4Optim.x_sym — Function

x_sym(mod::Model)

Return variable model parameter names.

Default

Returns mod.x_sym::Vector{Symbol}.

VP4Optim.x — Function

x(mod::Model)

Returns variable model parameters, according to the order specified in mod.x_ind.

Default

Returns mod.val[mod.x_ind]::Vector{Float64}.

VP4Optim.x! — Method

x!(mod::Model, x_syms::AbstractArray, x_vals::AbstractArray)

Resets those variable parameters, which are specified in x_syms, by the values in x_vals.

Default

Copies the values from x_vals into the associated locations.
Subsequently calls x_changed! to trigger optional secondary actions.
Returns nothing.

VP4Optim.x! — Method

x!(mod::Model, new_x::AbstractArray)

Resets all variable parameters by the values in new_x.

Default

Copies the values from new_x into val[x_ind] (in this order).
Subsequently calls x_changed! to trigger optional secondary actions.
Returns nothing.

VP4Optim.par_sym — Function

par_sym(mod::Model)

Returns fixed model parameter names.

Default

Returns mod.par_sym::Vector{Symbol}.

VP4Optim.par — Function

par(mod::Model)

Returns fixed model parameters, according to the order specified in mod.par_ind.

Default

Returns mod.val[mod.par_ind]::Vector{Float64}.

VP4Optim.par! — Method

par!(mod::Model, p_syms::AbstractArray, p_vals::AbstractArray)

Resets those fixed parameters, which are specified in par_syms, by the values in p_vals.

Default

Copies the values from p_vals into the associated locations.
Subsequently calls par_changed! to trigger optional secondary actions.
Returns nothing.

VP4Optim.par! — Method

par!(mod::Model, new_par::AbstractArray)

Resets all fixed parameters by the values in new_par.

Default

Copies the values from new_par into val[par_ind] (in this order).
Subsequently calls par_changed! to trigger optional secondary actions.
Returns nothing.

VP4Optim.y — Function

y(mod::Model)

Returns the actual data vector.

Default

Returns mod.y::SVector{Ny, T}.

VP4Optim.y! — Function

y!(mod::Model{Ny,Nx,Nc,T}, new_y::AbstractArray) where {Ny,Nx,Nc,T}

Sets new data values.

Default

Resets mod.y::SVector{Ny, T} with the content of new_y.
Calculates the squared magnitude of mod.y and stores the result in mod.y2::Float64.
Returns nothing.

Remarks

In most cases, it is safer to call set_data! instead of y!.

Note

Changing the data does not require to generate a new model instance!
To be on the safe side, use set_data! instead of y!.

VP4Optim.set_data! — Function

set_data!(mod::Model, new_data::AbstractArray)

Sets new data values safely.

Default

Simply calls y!.

Remarks

Can be used to perform preprocessing of data, if required by specific models.
This is the recommended method to (re)set new data.

VARPRO routines

These methods provide an interface to VARPRO, as described in Variable Projection.

VP4Optim.A — Function

A(::Model)

Return VARPRO matrix A.

Default

Returns mod.A, if the field exists and throws an error otherwise.

Remarks

Proper functioning of the method is mandatory. If the field mod.A does not exist, the method has to be implemented for the specific model.
Can be replaced by model-specific implementation, if the field mod.A does not exist.
In that scenario, the methods x_changed! and par_changed! trigger updates of mod.A.
if the models exhibit redundancy , the return type can be

VP4Optim.B — Function

B(mod::Model)

Returns matrix B::SMatrix{Nc,Nc}

Remarks

Returns A(mod)' * A(mod) by default.

VP4Optim.b — Function

b(mod::Model)

Returns vector b::SVector{Nc}

Remarks

Returns A(mod)' * y(mod) by default.

VP4Optim.Bb! — Function

Bb!(mod::Model)

Return matrix B = A' * A and vector b = A' * y.

Default

Direct calculation, based upon methods A and y.

Remarks

Mandatory routine for the calculation of χ² (and its derivatives)
Can be replaced by model-specific implementation, to improve the performance.
Expected to return (B, b)::Tuple.
For general models, the return types could be B::SMatrix{Nc,Nc,T} and b::SVector{Nc,T}
In most cases, model-specific implementations will be more efficient though.

Note

As long as y^2 - b' * (B \ b) gives the correct result (χ2), everything should be fine.
The same is true for the partial derivatives.

In case of doubt, look at how f, fg! and fgh! are implemented in the source code.

VP4Optim.c — Function

c(mod::Model)

Return VARPRO vector c.

Default

Gets B and b from Bb! and calculates generic solution c = B \ b.

Remarks

Can be replaced by model-specific implementation, if desired (e.g. for performance improvements).

VP4Optim.y_model — Function

y_model(mod::Model)

Compute model prediction A(x) * c.

Default

Calculates the product of the methods A and c

Remarks

Return type == SVector{Ny,T}
Can be used to check the model or generate synthetic data.
If necessary, a model-specific implementation may be needed. (e.g., if the method A has not be implemented.)

VP4Optim.χ2 — Function

χ2(mod::Model)

Return χ² = y² - b' * (B \ b) of actual model.

Default

Uses mod.y2 and (B, b) from Bb! to directly calculate the expression.

Interface for Optim.jl

Numerical minimization can be facilitated by the use of powerful optimization libraries, such as Optim.jl. To this end, VP4Optim provides some convenience functions, which provide interfaces to $χ²$ and its partial derivatives (up to second order) in a form as typically expected by the optimization libraries. Also a Hessian-based preconditioner is available.

VP4Optim.f — Method

f(mod::Model)

Return function f of argument x to be minimized, as expected by Optim.jl

Remark

Returns anonymous function x -> ... (cf. Optim.jl)
Depends on Bb!.

VP4Optim.fg! — Method

fg!(mod::Model)

Return function fg! of three arguments (F, G, x) as expected by Optim.jl.

Remark

Returns anonymous function (F, G, x) -> ... (cf. Optim.jl)
Depends on Bb! and ∂Bb!.

VP4Optim.fgh! — Method

fgh!(mod::Model)

Return function fgh! of four arguments (F, G, H, x) as expected by Optim.jl.

Remark

Returns anonymous function (F, G, H, x) -> ... (cf. Optim.jl)
Depends on Bb!, ∂Bb! and ∂∂Bb!.

VP4Optim.P — Function

P(mod::Model{Ny,Nx,Nc,T}, x) where {Ny,Nx,Nc,T}

Returns Hessian of $χ²(x)$.

Remark

Can be used as preconditioner, as expected by Optim.jl.

Model specific

VP4Optim.x_changed! — Function

x_changed!(::Model)

Informs user-defined model that x has changed.

Default

Does nothing.

Remarks

Can be used to recalculate any auxiliary model variable (such as A), which depends on x.

VP4Optim.par_changed! — Function

par_changed!(::Model)

Informs user-defined model that par has changed.

Default

Does nothing.

Remarks

Can be used to recalculate any auxiliary model variable (such as A), which depends on par.

VP4Optim.∂Bb! — Function

∂Bb!(::Model)

Returns up to first order partial derivatives with respect to x.

Default

None, must be supplied by the user.

Remarks

Required for first and second order optimization techniques.
Returns (B, b, ∂B, ∂b)::Tuple
(B, b) as returned from Bb!
Types: ∂B::SVector{Nx, SMatrix{Nc,Nc,T}} and ∂b::SVector{Nx, SVector{Nc,T}}
(or anything, which works with the implementation of fg! and fgh!)

VP4Optim.∂∂Bb! — Function

∂∂Bb!(::Model)

Returns up to second order partial derivatives with respect to x

Default

None, must be supplied by the user.

Remarks

Required for second order optimization techniques.
Returns (B, b, ∂B, ∂b, ∂∂B, ∂∂b)::Tuple
(B, b, ∂B, ∂b) as returned from ∂Bb!
Types: ∂∂B::SMatrix{Nx, Nx, SMatrix{Nc,Nc,T}} and ∂∂b::SMatrix{Nx, Nx, SVector{Nc,T}}
(or anything, which works with the implementation of fgh!)

Model testing

VP4Optim.check_model — Function

check_model(pars, vals, c_, data_;
what=(:consistency, :derivatives, :optimization),
small=sqrt(eps()),
x0=[], lx=[], ux=[], x_scale=[],
precon=true,
visual=false,
rng=MersenneTwister(),
Hessian=true,
log10_rng=range(-6, -3, 10),
min_slope=0.9)

Tests, which any specific model should pass.

Arguments

pars::ModPar{T}: Model parameters for constructor create_model(pars).
vals::Vector{Float64}: All nonlinear parameters, the model depends on. As defined in the Model field val.
c_::Vector{Nc, T}: Linear cofficients, the model depends on.
data_::AbstractArray: Data, corresponding to the true parameters vals and c_, in a format suitable for set_data!.
what::Tuple{Symbol}: Tests to be performed (see below).
small::Float64: Required as accuracy criterion.
x0::Vector{Float64}: Starting point for optimization and location, where derivatives are tested.
lx::Vector{Float64}: Lower bound of optimization
ux::Vector{Float64}: Upper bound of optimization
x_scale::Vector{Float64}: Scaling vector, such that δx = randn(size(x)) .* x_scale becomes reasonable
precon::Bool: Test optimization with and without preconditioner.
visual::Bool: If true also generate double-logarithmic plots for the derivative tests.
rng::MersenneTwister: Allows to pass a unique seed (e.g. MersenneTwister(42)) for reproducible testing.
Hessian::Bool: Should be set to false, if the model does not implement second order derivatives.
log10_rng::AbstractVector: logarithmic range for derivative testing
min_slope::Float64: minimal derivative slope on log-log plot

Remark

Tests are performed for every x_sym ⊆ sym.
Returns a dictionary with detailed information about the test results.
:consistency ∈ what: Several basic tests (parameters, names, correct model values)
:derivatives ∈ what: Check first and second order partial derivatives at x0.
:optimization ∈ what: Minimize model with x0 as starting point and bounds lx and ux.
An example application can be found in test_BiExpDecay.jl. This should also work as a template, how to perform tests on own models.

1In Variable Projection, $\bm{x}$ referred to this variable part only, while any fixed parameters were absorbed in the actual definition of the matrix $\bm{A}$.
2Here, we assume a fixed data type for the specific model.