Add WyNDA online estimator (#631)

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>

Author: ChrisRackauckas-Claude

docs/src/libs/datadrivensparse/sparse_regression.md

@@ -52,6 +52,7 @@ Where the matrix of evaluated basis elements $\varPhi_y \in \mathbb R^{\lvert \v
 STLSQ
 ADMM
 SR3
+WyNDA
 ImplicitOptimizer
 ```
 

lib/DataDrivenSparse/src/DataDrivenSparse.jl

@@ -123,6 +123,9 @@ export ADMM
 include("algorithms/SR3.jl")
 export SR3
 
+include("algorithms/WyNDA.jl")
+export WyNDA
+
 include("algorithms/Implicit.jl")
 export ImplicitOptimizer
 

lib/DataDrivenSparse/src/algorithms/WyNDA.jl

@@ -0,0 +1,102 @@
+"""
+$(TYPEDEF)
+
+`WyNDA` implements the Wide-Array of Nonlinear Dynamics Approximation update as an
+online recursive least-squares estimator with exponential forgetting. Given a
+matrix of evaluated basis functions `X` and targets `Y`, it updates the
+coefficient matrix sample-by-sample so that `Y ≈ coefficients * X`.
+
+The forgetting factor `λ` controls how quickly older samples are discounted.
+Values close to one recover a batch least-squares-like fit, while smaller values
+adapt faster to parameter drift.
+
+# Fields
+
+$(FIELDS)
+
+# Example
+
+```julia
+opt = WyNDA()
+opt = WyNDA(0.998)
+opt = WyNDA(0.998; initial_covariance = 1.0e4)
+```
+"""
+struct WyNDA{T <: Number, C, IC} <: AbstractSparseRegressionAlgorithm
+    """Exponential forgetting factor."""
+    λ::T
+    """Initial inverse covariance scale or matrix."""
+    initial_covariance::C
+    """Optional initial coefficient vector or matrix."""
+    initial_coefficients::IC
+
+    function WyNDA(
+            λ::T = 1.0;
+            initial_covariance::C = 1.0e6,
+            initial_coefficients::IC = nothing
+        ) where {T <: Number, C, IC}
+        @assert zero(T) < λ <= one(T) "Forgetting factor λ must satisfy 0 < λ <= 1"
+        if initial_covariance isa Number
+            @assert initial_covariance > zero(initial_covariance) "Initial covariance must be positive"
+        end
+        return new{T, C, IC}(λ, initial_covariance, initial_coefficients)
+    end
+end
+
+Base.summary(::WyNDA) = "WyNDA"
+
+function _initial_covariance(alg::WyNDA, ::Type{T}, n_features::Int) where {T}
+    covariance = alg.initial_covariance
+    if covariance isa Number
+        return Matrix{T}(I, n_features, n_features) .* T(covariance)
+    end
+    @assert size(covariance) == (n_features, n_features) "Initial covariance matrix size must match the number of features"
+    return Matrix{T}(covariance)
+end
+
+function _initial_coefficients(
+        alg::WyNDA, ::Type{T}, n_targets::Int, n_features::Int
+    ) where {T}
+    coefficients = alg.initial_coefficients
+    if coefficients === nothing
+        return zeros(T, n_targets, n_features)
+    end
+    if coefficients isa AbstractVector
+        @assert n_targets == 1 "Initial coefficient vectors are only valid for single-target problems"
+        @assert length(coefficients) == n_features "Initial coefficient vector length must match the number of features"
+        return reshape(T.(coefficients), 1, n_features)
+    end
+    @assert size(coefficients) == (n_targets, n_features) "Initial coefficient matrix size must match targets by features"
+    return Matrix{T}(coefficients)
+end
+
+function (alg::WyNDA)(
+        X::AbstractMatrix, Y::AbstractVecOrMat;
+        options::DataDrivenCommonOptions = DataDrivenCommonOptions(),
+        kwargs...
+    )
+    Y_matrix = Y isa AbstractVector ? reshape(Y, 1, :) : Y
+    @assert size(X, 2) == size(Y_matrix, 2) "X and Y must have the same number of observations"
+
+    n_features, n_observations = size(X)
+    n_targets = size(Y_matrix, 1)
+    T = float(promote_type(eltype(X), eltype(Y_matrix), typeof(alg.λ)))
+
+    coefficients = _initial_coefficients(alg, T, n_targets, n_features)
+    covariance = _initial_covariance(alg, T, n_features)
+    λ = T(alg.λ)
+
+    for k in 1:n_observations
+        basis_values = view(X, :, k)
+        targets = view(Y_matrix, :, k)
+        covariance_basis = covariance * basis_values
+        denominator = λ + dot(basis_values, covariance_basis)
+        gain = covariance_basis ./ denominator
+        residual = targets .- coefficients * basis_values
+        coefficients .+= residual * adjoint(gain)
+        covariance .-= gain * adjoint(covariance_basis)
+        covariance ./= λ
+    end
+
+    return coefficients, alg.λ, n_observations
+end

lib/DataDrivenSparse/test/Core/sparse_linear_solve.jl

@@ -97,6 +97,34 @@ end
     end
 end
 
+@testset "WyNDA online estimator" begin
+    rng = StableRNG(20260708)
+    t = 0.0:0.02:8.0
+    X = permutedims(
+        reduce(
+            hcat,
+            (
+                ones(length(t)), sin.(0.7 .* t), cos.(1.3 .* t),
+                exp.(-0.2 .* t), randn(rng, length(t)),
+            )
+        )
+    )
+    A = [0.5 -1.2 0.0 0.75 0.0; -0.25 0.0 1.5 0.0 0.1]
+    Y = A * X
+
+    coefficients, λ, iterations = WyNDA(1.0; initial_covariance = 1.0e8)(X, Y)
+
+    @test λ == 1.0
+    @test iterations == size(X, 2)
+    # Finite initial covariance is equivalent to a small ridge prior in RLS.
+    @test coefficients ≈ A atol = 1.0e-6
+
+    problem = DirectDataDrivenProblem(X, Y)
+    result = @test_nowarn solve(problem, WyNDA(1.0; initial_covariance = 1.0e8))
+    @test result isa DataDrivenSolution
+    @test result.residuals <= 1.0e-8
+end
+
 # Issue #564: Test that solve doesn't throw MethodError when coefficients are all zero
 # This can happen with very small data values or high regularization
 @testset "Zero coefficients handling (Issue #564)" begin