From 6546324d887c938a10abd860624241ec3eee190a Mon Sep 17 00:00:00 2001 From: MartinuzziFrancesco Date: Tue, 30 Jun 2026 13:37:54 +0200 Subject: [PATCH 1/2] Add conceptors (Jaeger 2014) Implement conceptors for recurrent neural networks after Jaeger (2014), "Controlling Recurrent Neural Networks by Conceptors" (arXiv:1403.3369), as a single src/conceptors.jl included in the module: - conceptor matrices C = R (R + a^-2 I)^-1, correlation, singular values, quota - aperture adaptation phi(C, gamma), reaperture, attenuation, optimal_aperture - Boolean algebra: conceptor_not / and / or (singular-safe), with unexported infix sugar not/and/or - Conceptor reservoir-computer wrapper + named conceptor library - load! (input-internalizing recurrent weights), generate (autonomous rollout), morph_conceptor (linear conceptor mixtures) - conceptor-filtered supervised training (store_conceptors!, train! method) Adds test/test_conceptors.jl (auto-discovered Core test, 42 tests), a docs example reproducing the morphing square (Fig. 2), an API page, and a bib entry. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K8VYJcvx5GAAd16NEjL8Jm --- docs/pages.jl | 2 + docs/src/api/conceptors.md | 58 ++ docs/src/examples/conceptors_morphing.md | 127 +++++ docs/src/refs.bib | 9 + src/ReservoirComputing.jl | 15 +- src/conceptors.jl | 643 +++++++++++++++++++++++ test/test_conceptors.jl | 131 +++++ 7 files changed, 983 insertions(+), 2 deletions(-) create mode 100644 docs/src/api/conceptors.md create mode 100644 docs/src/examples/conceptors_morphing.md create mode 100644 src/conceptors.jl create mode 100644 test/test_conceptors.jl diff --git a/docs/pages.jl b/docs/pages.jl index a15e1791..ee608e18 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -14,10 +14,12 @@ pages = [ ], "Examples" => Any[ "Building a model to add to ReservoirComputing.jl" => "examples/model_es2n.md", + "Morphing patterns with conceptors" => "examples/conceptors_morphing.md", ], "API Documentation" => Any[ "Layers" => "api/layers.md", "Models" => "api/models.md", + "Conceptors" => "api/conceptors.md", "Utilities" => "api/utils.md", "Train" => "api/train.md", "Predict" => "api/predict.md", diff --git a/docs/src/api/conceptors.md b/docs/src/api/conceptors.md new file mode 100644 index 00000000..79d8012c --- /dev/null +++ b/docs/src/api/conceptors.md @@ -0,0 +1,58 @@ +# Conceptors + +Conceptors after Jaeger (2014), [Jaeger2014conceptors](@cite). See the +[Morphing patterns with conceptors](@ref) example for an end-to-end walkthrough. + +## Conceptor matrices + +```@docs + conceptor_matrix + correlation_matrix + conceptor_from_states + conceptor_singular_values + quota +``` + +## Aperture adaptation + +```@docs + aperture_adapt + adapt_singular_value + reaperture + attenuation + optimal_aperture +``` + +## Boolean algebra + +```@docs + conceptor_not + conceptor_and + conceptor_or +``` + +## The `Conceptor` wrapper and its library + +```@docs + Conceptor + has_conceptor + get_conceptor + store_conceptor! + set_active_conceptor + active_conceptor +``` + +## Loading, generation, and morphing + +```@docs + load! + generate + morph_conceptor + ridge_map +``` + +## Conceptor-filtered training + +```@docs + store_conceptors! +``` diff --git a/docs/src/examples/conceptors_morphing.md b/docs/src/examples/conceptors_morphing.md new file mode 100644 index 00000000..2a325174 --- /dev/null +++ b/docs/src/examples/conceptors_morphing.md @@ -0,0 +1,127 @@ +# Morphing patterns with conceptors + +This example reproduces the morphing square from Jaeger (2014), *Controlling +Recurrent Neural Networks by Conceptors* ([Jaeger2014conceptors](@cite), +arXiv:1403.3369, Figure 2). A single reservoir is *loaded* with four patterns; +their [`conceptor`](@ref conceptor_matrix)s are then linearly **morphed** and used +to drive the reservoir autonomously, generating a continuum of patterns that +interpolate between, and extrapolate beyond, the four originals. + +A conceptor is a positive semidefinite matrix ``C = R (R + \alpha^{-2} I)^{-1}`` +built from the correlation matrix ``R`` of a reservoir's driven states. Inserted +into the autonomous update ``x(n) = C\,\tanh(W x(n-1) + b)``, it constrains the +reservoir dynamics to the subspace excited by one loaded pattern, so the reservoir +re-generates that pattern. A linear combination of conceptors morphs between them. + +## The four driving patterns + +Following the original demonstration we use two sines of slightly different +period and two minor variations of a 5-periodic random pattern, all in +``[-0.9, 0.9]``. The two period-5 vectors are the ones from the published figure. + +```@example morphing +using ReservoirComputing +using Random +using Plots + +const PERIOD_1 = 8.8342522 +const PERIOD_2 = 9.8342522 +const DRIVER_LEN = 1500 + +# period-5 patterns, sampled as rp[mod(n,5)+1] (values from Jaeger 2014, Fig. 1B) +const RAND5 = [0.16, 0.9, -0.9, -0.21, -0.55] +const PERTURB5 = [0.10, 0.9, -0.9, -0.65, -0.54] +period5(rp, n) = rp[mod(n, 5) + 1] + +n = 1:DRIVER_LEN +patterns = [ + :s1 => reshape(sin.(2π .* n ./ PERIOD_1), 1, :), + :s2 => reshape(sin.(2π .* n ./ PERIOD_2), 1, :), + :r3 => reshape([period5(RAND5, i) for i in n], 1, :), + :r4 => reshape([period5(PERTURB5, i) for i in n], 1, :), +] +nothing # hide +``` + +## Building and loading the reservoir + +We build a 100-unit [`ESN`](@ref) with the scalings used in the report +(spectral radius 1.5, input scaling 1.5, bias scaling 0.2) and wrap it in a +[`Conceptor`](@ref). [`load!`](@ref) drives the reservoir with each pattern, +stores a conceptor for it, and recomputes the recurrent weights into an +input-internalizing matrix so the reservoir can run autonomously. + +```@example morphing +rng = Xoshiro(3) + +bias_init(r, dims...) = 0.2f0 .* randn(r, Float32, dims...) +input_init(r, dims...) = 1.5f0 .* randn(r, Float32, dims...) +res_init(r, dims...) = rand_sparse(r, dims...; radius = 1.5f0, sparsity = 0.1f0) + +esn = ESN(1, 100, 1; use_bias = true, init_bias = bias_init, + init_input = input_init, init_reservoir = res_init) +concept = Conceptor(esn) +ps = initialparameters(rng, concept) +st = initialstates(rng, concept) + +ps, st = load!(rng, concept, patterns, ps, st; aperture = 4.0, washout = 500) +nothing # hide +``` + +## Morphing across the square + +For a grid of mixing coordinates ``a, b \in \{-0.5, \dots, 1.5\}`` we form the +morphed conceptor ``M = \mu_1 C_{s1} + \mu_2 C_{s2} + \mu_3 C_{r3} + \mu_4 C_{r4}`` +with +```math +\mu_1 = (1-a)\,b,\quad \mu_2 = a\,b,\quad \mu_3 = (1-a)(1-b),\quad \mu_4 = a\,(1-b), +``` +(the coefficients always sum to one), and run the reservoir autonomously under +``M`` with [`generate`](@ref). All panels share one start state so their phases are +comparable. + +```@example morphing +morph_weights(a, b) = (; s1 = (1 - a) * b, s2 = a * b, + r3 = (1 - a) * (1 - b), r4 = a * (1 - b)) + +grid = collect(-0.5:0.25:1.5) +x0 = rand(Xoshiro(4), 100) +prototypes = Set([(0.0, 1.0), (1.0, 1.0), (0.0, 0.0), (1.0, 0.0)]) + +plt = plot(layout = (length(grid), length(grid)), size = (1000, 1000), + legend = false, framestyle = :box, ticks = false, link = :all) + +for (ib, b) in enumerate(grid), (ia, a) in enumerate(grid) + M = morph_conceptor(st, morph_weights(a, b)) + Y, _ = generate(concept, ps, st; conceptor = M, steps = 15, + washout = 190, init_state = x0) + k = (ib - 1) * length(grid) + ia # row ib (top = b = -0.5), col ia + proto = (a, b) in prototypes + plot!(plt, subplot = k, vec(Y); color = :black, linewidth = proto ? 2.5 : 1.2, + ylims = (-1, 1), background_color_subplot = proto ? :gray88 : :white) +end +plt +``` + +The four shaded panels are the loaded prototypes: the two sines sit on the +``b = 1`` row, the two period-5 patterns on the ``b = 0`` row. The inner block +(``a, b \in [0, 1]``) interpolates smoothly between them; the outer panels +extrapolate, producing distorted, higher-frequency variants. + +## Single-pattern recall and the conceptor algebra + +Loading also lets the reservoir recall a single pattern by name, and the stored +conceptors support an aperture adaptation ([`aperture_adapt`](@ref)) and a Boolean +algebra ([`conceptor_and`](@ref), [`conceptor_or`](@ref), [`conceptor_not`](@ref)). + +```@example morphing +# autonomous recall of the first sine +Yrec, _ = generate(concept, ps, st; conceptor = :s1, steps = 100, washout = 200) + +# conceptor algebra: shared vs. combined subspaces of the two sines +Cs1, Cs2 = get_conceptor(st, :s1), get_conceptor(st, :s2) +(quota(conceptor_and(Cs1, Cs2)), quota(Cs1), quota(conceptor_or(Cs1, Cs2))) +``` + +The triple is increasing — `AND` keeps only the shared directions, `OR` spans the +union — illustrating the lattice structure of the conceptor algebra. diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 0f3946ab..2882e3e3 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -1,3 +1,12 @@ +@article{Jaeger2014conceptors, + title = {Controlling Recurrent Neural Networks by Conceptors}, + url = {https://arxiv.org/abs/1403.3369}, + DOI = {10.48550/arXiv.1403.3369}, + journal = {arXiv preprint arXiv:1403.3369}, + author = {Jaeger, Herbert}, + year = {2014}, +} + @article{Lu2017, title = {Reservoir observers: Model-free inference of unmeasured variables in chaotic systems}, volume = {27}, diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl index 3bec529a..f38340a4 100644 --- a/src/ReservoirComputing.jl +++ b/src/ReservoirComputing.jl @@ -2,12 +2,13 @@ module ReservoirComputing using ArrayInterface: ArrayInterface using ConcreteStructs: @concrete -using LinearAlgebra: eigvals, I, qr, Diagonal, diag, mul!, Symmetric, norm +using LinearAlgebra: eigvals, eigen, I, qr, Diagonal, diag, mul!, Symmetric, norm, + svd, svdvals, nullspace, pinv, tr using LuxCore: AbstractLuxLayer, AbstractLuxContainerLayer, AbstractLuxWrapperLayer, setup, apply, replicate import LuxCore: initialparameters, initialstates, statelength, outputsize using NNlib: tanh_fast -using Random: Random, AbstractRNG, randperm +using Random: Random, AbstractRNG, randperm, randn using Static: StaticBool, StaticSymbol, True, False, static, known, StaticInteger using Reexport: Reexport, @reexport using WeightInitializers: WeightInitializers, DeviceAgnostic, PartialFunction, Utils, @@ -59,6 +60,8 @@ include("models/lifesn.jl") include("models/ngrc.jl") include("models/rmnesn.jl") include("models/rmnresesn.jl") +#conceptors +include("conceptors.jl") #extensions include("extensions/reca.jl") @@ -85,6 +88,14 @@ export polynomial_monomials, chebyshev_monomials, predict, QRSolver, resetcarry! export AdditiveEIESN, DeepESN, DelayESN, EIESN, ES2N, ESN, EuSN, HybridESN, InputDelayESN, LIFESN, ResESN, StateDelayESN, SVESM export NGRC export RMNESN, RMNResESN +#conceptors +export Conceptor, correlation_matrix, conceptor_matrix, conceptor_from_states, + conceptor_singular_values, quota +export aperture_adapt, adapt_singular_value, reaperture, attenuation, optimal_aperture +export conceptor_not, conceptor_and, conceptor_or +export has_conceptor, get_conceptor, store_conceptor!, store_conceptors!, + set_active_conceptor, active_conceptor +export ridge_map, load!, generate, morph_conceptor #ext export RECACell, RECA export RandomMapping, RandomMaps diff --git a/src/conceptors.jl b/src/conceptors.jl new file mode 100644 index 00000000..1c5edc17 --- /dev/null +++ b/src/conceptors.jl @@ -0,0 +1,643 @@ +# Conceptors for recurrent neural networks, after Jaeger (2014), "Controlling +# Recurrent Neural Networks by Conceptors" (arXiv:1403.3369). +# +# A conceptor is a positive semidefinite matrix C = R (R + α^{-2} I)^{-1} derived +# from the correlation matrix R of a reservoir's driven states. It acts as a soft +# projector onto the linear subspace a driving pattern excites in the reservoir. +# Conceptors can be aperture-adapted, combined with a Boolean algebra, loaded into +# a reservoir to autonomously regenerate the loaded patterns, and morphed to +# interpolate or extrapolate between patterns. +# +# Conceptor matrices are always computed and stored in `Float64`: the defining +# operation is a matrix inversion whose conditioning degrades quickly for the +# near-singular correlation matrices produced by periodic drivers, so conceptor +# precision is decoupled from the (often `Float32`) reservoir precision. + +const ConceptorMatrix = Matrix{Float64} + +# ====================================================================== +# Core conceptor matrices +# ====================================================================== + +@doc raw""" + correlation_matrix(states) -> Matrix{Float64} + +State correlation matrix ``R = X X^\top / L`` for a reservoir state collection +`states` of size ``N \times L`` (one state per column). This is the matrix whose +regularized-identity map defines a conceptor. +""" +function correlation_matrix(states::AbstractMatrix{<:Real}) + X = Float64.(states) + return (X * X') / size(X, 2) +end + +@doc raw""" + conceptor_matrix(R, aperture) -> Matrix{Float64} + +Conceptor ``C = R (R + \alpha^{-2} I)^{-1}`` derived from a correlation matrix `R` +with aperture ``\alpha`` = `aperture` (Jaeger 2014, Eq. 7). The result is +symmetric positive semidefinite with singular values in ``[0, 1)``. +""" +function conceptor_matrix(R::AbstractMatrix{<:Real}, aperture::Real) + aperture > 0 || throw(ArgumentError("aperture must be positive, got $aperture")) + Rd = Float64.(R) + inv_sq = 1.0 / aperture^2 + return Rd / (Rd + inv_sq * I) +end + +""" + conceptor_from_states(states, aperture) -> Matrix{Float64} + +Convenience composition of [`correlation_matrix`](@ref) and +[`conceptor_matrix`](@ref): the conceptor characterizing a reservoir state cloud +`states` (size `N × L`) at the given `aperture`. +""" +function conceptor_from_states(states::AbstractMatrix{<:Real}, aperture::Real) + return conceptor_matrix(correlation_matrix(states), aperture) +end + +""" + conceptor_singular_values(C) -> Vector{Float64} + +Singular values of a (symmetric) conceptor matrix `C`, in descending order. For a +conceptor these coincide with its eigenvalues and lie in `[0, 1]`. +""" +function conceptor_singular_values(C::AbstractMatrix{<:Real}) + return svdvals(Float64.(C)) +end + +""" + quota(C) -> Float64 + +Mean singular value of a conceptor, `tr(C) / N`. Jaeger's "quota" measures the +fraction of reservoir state space the conceptor leaves open (0 = point, 1 = all). +""" +function quota(C::AbstractMatrix{<:Real}) + return tr(Float64.(C)) / size(C, 1) +end + +# ====================================================================== +# Aperture adaptation and aperture selection +# ====================================================================== + +@doc raw""" + adapt_singular_value(s, γ) -> Float64 + +Aperture-adapted singular value ``s_\gamma`` for an input singular value +``s \in [0, 1]`` and adaptation factor ``\gamma \in [0, \infty]`` (Jaeger 2014, +Prop. 3, Eq. 19). Hard values `s = 0` and `s = 1` are fixed points; intermediate +values are pushed toward 0 as ``\gamma \to 0`` and toward 1 as ``\gamma \to \infty``. +""" +function adapt_singular_value(s::Float64, γ::Float64) + (s <= 0.0) && return 0.0 + (s >= 1.0) && return 1.0 + iszero(γ) && return 0.0 + isinf(γ) && return 1.0 + return s / (s + (1.0 - s) / γ^2) +end + +@doc raw""" + aperture_adapt(C, γ) -> Matrix{Float64} + +Aperture adaptation ``\varphi(C, \gamma)`` of a conceptor `C` by factor +``\gamma \in [0, \infty]`` (Jaeger 2014, Def. 3). Equivalent to multiplying the +aperture of `C` by ``\gamma``; in particular ``C(R, \alpha) = \varphi(C(R, 1), \alpha)``. +Applied via the eigendecomposition so the limiting cases ``\gamma = 0`` (hard +zero) and ``\gamma = \infty`` (hardening to a projector) are exact. +""" +function aperture_adapt(C::AbstractMatrix{<:Real}, γ::Real) + γ >= 0 || throw(ArgumentError("aperture factor γ must be ≥ 0, got $γ")) + F = eigen(Symmetric(Float64.(C))) + s_adapted = adapt_singular_value.(clamp.(F.values, 0.0, 1.0), Float64(γ)) + return F.vectors * Diagonal(s_adapted) * F.vectors' +end + +""" + reaperture(C, from_aperture, to_aperture) -> Matrix{Float64} + +Re-express a conceptor `C` that currently has aperture `from_aperture` so that it +has aperture `to_aperture`, using `φ(C, to/from)` (Jaeger 2014, Prop. 5). +""" +function reaperture(C::AbstractMatrix{<:Real}, from_aperture::Real, to_aperture::Real) + from_aperture > 0 || throw(ArgumentError("from_aperture must be positive")) + return aperture_adapt(C, to_aperture / from_aperture) +end + +@doc raw""" + attenuation(C, W, b; steps, washout, init_state, rng) -> Float64 + +Attenuation ``a_C = E[\|z(n) - x(n)\|^2] / E[\|z(n)\|^2]`` of a loaded reservoir +(recurrent weights `W`, bias `b`) run autonomously under conceptor `C` +(Jaeger 2014, Eq. 23), where ``z(n) = \tanh(W x(n-1) + b)`` is the unconstrained +update and ``x(n) = C z(n)`` is the conceptor-constrained state. The attenuation +is the fraction of reservoir signal energy suppressed by `C`; as a function of +aperture it passes through a minimum at the best-reconstructing aperture. +""" +function attenuation( + C::AbstractMatrix{<:Real}, W::AbstractMatrix{<:Real}, b::AbstractVector{<:Real}; + steps::Int = 500, washout::Int = 200, + init_state::Union{Nothing, AbstractVector} = nothing, + rng::AbstractRNG = Random.default_rng() + ) + Cd = Float64.(C) + Wd = Float64.(W) + bd = Float64.(b) + N = size(Wd, 1) + x = init_state === nothing ? 0.5 .* randn(rng, N) : Float64.(init_state) + + num = 0.0 + den = 0.0 + for n in 1:(washout + steps) + z = tanh.(Wd * x .+ bd) + x = Cd * z + if n > washout + num += sum(abs2, z .- x) + den += sum(abs2, z) + end + end + return num / den +end + +""" + optimal_aperture(R, W, b, apertures; kwargs...) -> (best_aperture, attenuations) + +Select the aperture minimizing the [`attenuation`](@ref) criterion over a grid +`apertures`. For each candidate `α` the conceptor `C(R, α)` is formed and the +loaded reservoir (`W`, `b`) is run autonomously to measure its attenuation. Returns +the minimizing aperture together with the full vector of attenuations (aligned with +`apertures`) so the characteristic trough can be inspected. `kwargs` are forwarded +to [`attenuation`](@ref). +""" +function optimal_aperture( + R::AbstractMatrix{<:Real}, W::AbstractMatrix{<:Real}, b::AbstractVector{<:Real}, + apertures::AbstractVector{<:Real}; kwargs... + ) + atts = [attenuation(conceptor_matrix(R, α), W, b; kwargs...) for α in apertures] + return apertures[argmin(atts)], atts +end + +# ====================================================================== +# Boolean operations on conceptors (Jaeger 2014, Section 3.9) +# ====================================================================== + +@doc raw""" + conceptor_not(C) -> Matrix{Float64} + +Logical negation ``\neg C = I - C`` of a conceptor (Jaeger 2014, Eq. 28). +Exchanges the roles of the directions the conceptor admits and suppresses. +""" +function conceptor_not(C::AbstractMatrix{<:Real}) + Cd = Float64.(C) + return Matrix(I, size(Cd)) - Cd +end + +# Orthonormal basis of the range of a symmetric PSD matrix, using a tolerance on +# the singular values to decide the numerical rank. +function _range_basis(C::Matrix{Float64}; tol::Float64 = 1.0e-12) + F = svd(C) + rank = count(>(tol), F.S) + return F.U[:, 1:rank] +end + +@doc raw""" + conceptor_and(C, B) -> Matrix{Float64} + +Logical conjunction ``C \wedge B`` of two conceptors (Jaeger 2014, Eq. 32). For +full-rank conceptors this equals ``(C^{-1} + B^{-1} - I)^{-1}``; the implementation +uses the range-intersection projector form so that singular conceptors are handled +robustly. The result admits exactly the reservoir directions admitted by *both* +`C` and `B`. +""" +function conceptor_and(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) + Cd = Float64.(C) + Bd = Float64.(B) + size(Cd) == size(Bd) || throw(DimensionMismatch("conceptors must have equal size")) + + # Basis of R(C) ∩ R(B) = N(P_{N(C)} + P_{N(B)}), via the null-space projectors. + UC0 = nullspace(Cd) + UB0 = nullspace(Bd) + M = UC0 * UC0' + UB0 * UB0' + Fm = svd(M) + rank_M = count(>(1.0e-12), Fm.S) + Bgk = Fm.U[:, (rank_M + 1):end] # basis of the range intersection + + if size(Bgk, 2) == 0 + return zeros(size(Cd)) # disjoint supports → all-zero conceptor + end + core = Bgk' * (pinv(Cd) + pinv(Bd) - I) * Bgk + return Bgk * (core \ Matrix(I, size(core))) * Bgk' +end + +@doc raw""" + conceptor_or(C, B) -> Matrix{Float64} + +Logical disjunction ``C \vee B`` of two conceptors (Jaeger 2014, Eq. 31), defined +via De Morgan's law ``C \vee B = \neg(\neg C \wedge \neg B)``. The result admits +every reservoir direction admitted by `C` or by `B`. +""" +function conceptor_or(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) + return conceptor_not(conceptor_and(conceptor_not(C), conceptor_not(B))) +end + +# Infix sugar mirroring the paper's ¬ / ∧ / ∨ notation. `∧` and `∨` are Julia infix +# operators; `¬` is a prefix-style function (`¬(C)`). Not exported, to avoid +# polluting the package namespace; use the named functions or `ReservoirComputing.:∧`. +const ¬ = conceptor_not +∧(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) = conceptor_and(C, B) +∨(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) = conceptor_or(C, B) + +# ====================================================================== +# The Conceptor reservoir-computer wrapper and its parameter/state plumbing +# ====================================================================== + +@doc raw""" + Conceptor(model) + +Wrap a ReservoirComputing `model` (e.g. an [`ESN`](@ref)) so that conceptor +matrices can be derived from it, stored by name, combined with the conceptor +algebra, and used to constrain autonomous generation (Jaeger 2014). + +The wrapped `model` provides the reservoir update; the `Conceptor` adds a +conceptor library to the model state. Parameters and states are nested under the +`model` field, leaving room for conceptor bookkeeping at the top level. +""" +@concrete struct Conceptor <: AbstractReservoirComputer{(:model,)} + model +end + +function initialparameters(rng::AbstractRNG, concept::Conceptor) + return (; model = initialparameters(rng, concept.model), readout = nothing) +end + +function initialstates(rng::AbstractRNG, concept::Conceptor) + return (; + model = initialstates(rng, concept.model), + conceptors = Dict{Symbol, ConceptorMatrix}(), + apertures = Dict{Symbol, Float64}(), + active_conceptor = nothing, + ) +end + +""" + has_conceptor(st, name) -> Bool + +Whether a conceptor called `name` has been stored in the state `st`. +""" +has_conceptor(st::NamedTuple, name::Symbol) = haskey(st.conceptors, name) + +""" + get_conceptor(st, name) -> Matrix{Float64} or nothing + +The stored conceptor matrix called `name`, or `nothing` if it is absent. +""" +get_conceptor(st::NamedTuple, name::Symbol) = get(st.conceptors, name, nothing) + +""" + store_conceptor!(st, name, C, aperture) -> st + +Record conceptor matrix `C` (formed at `aperture`) under `name` in the state's +conceptor library. Mutates the library dictionaries in place and returns `st`. +""" +function store_conceptor!(st::NamedTuple, name::Symbol, C::AbstractMatrix{<:Real}, aperture::Real) + st.conceptors[name] = Float64.(C) + st.apertures[name] = Float64(aperture) + return st +end + +""" + set_active_conceptor(st, name) -> st + +Mark the conceptor `name` (or `nothing`) as the active one, returning an updated +state. The active conceptor is the default constraint for generation. +""" +function set_active_conceptor(st::NamedTuple, name::Union{Symbol, Nothing}) + name !== nothing && !has_conceptor(st, name) && throw(KeyError(name)) + return merge(st, (; active_conceptor = name)) +end + +""" + active_conceptor(st) -> Matrix{Float64} + +The currently active conceptor matrix. Throws if none has been set. +""" +function active_conceptor(st::NamedTuple) + name = st.active_conceptor + name === nothing && throw(ArgumentError("no active conceptor; call set_active_conceptor")) + return st.conceptors[name] +end + +# Resolve a conceptor argument that is either an explicit matrix or a stored name. +resolve_conceptor(st::NamedTuple, conceptor::AbstractMatrix{<:Real}) = Float64.(conceptor) +function resolve_conceptor(st::NamedTuple, conceptor::Symbol) + C = get_conceptor(st, conceptor) + C === nothing && throw(KeyError(conceptor)) + return C +end + +""" + collectstates(concept::Conceptor, signal, ps, st; conceptor=nothing) -> (states, st) + +Stream `signal` through the wrapped reservoir and collect its state sequence. When +`conceptor` (a stored `Symbol` or an explicit matrix) is given, the collected +states are filtered through it, `C x(n)`, realizing the conceptor as a state +modifier. +""" +function collectstates( + concept::Conceptor, signal::AbstractMatrix, ps::NamedTuple, st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix, Nothing} = nothing + ) + states, st_model = collectstates(concept.model, signal, ps.model, st.model) + new_st = merge(st, (; model = st_model)) + conceptor === nothing && return states, new_st + return resolve_conceptor(new_st, conceptor) * states, new_st +end + +""" + addreadout!(concept::Conceptor, W, ps, st) -> (ps, st) + +Install trained readout weights `W` into the wrapped model's readout layer +(`ps.model.readout.weight`), keeping the generation and classification paths on a +single readout convention. +""" +function addreadout!(concept::Conceptor, W::AbstractMatrix, ps::NamedTuple, st::NamedTuple) + return set_readout_weight(ps, W), st +end + +# ====================================================================== +# Parameter-tree access for the wrapped model (ps.model.{reservoir,readout}) +# ====================================================================== + +reservoir_params(ps::NamedTuple) = ps.model.reservoir +readout_params(ps::NamedTuple) = ps.model.readout + +function reservoir_bias(ps::NamedTuple) + rps = reservoir_params(ps) + return haskey(rps, :bias) ? Float64.(vec(rps.bias)) : zeros(Float64, size(rps.reservoir_matrix, 1)) +end + +function set_reservoir_matrix(ps::NamedTuple, W::AbstractMatrix) + reservoir = merge(ps.model.reservoir, (; reservoir_matrix = W)) + return merge(ps, (; model = merge(ps.model, (; reservoir = reservoir)))) +end + +function set_readout_weight(ps::NamedTuple, W::AbstractMatrix) + readout = merge(ps.model.readout, (; weight = W)) + return merge(ps, (; model = merge(ps.model, (; readout = readout)))) +end + +@doc raw""" + ridge_map(features, targets, reg) -> Matrix + +Tikhonov-regularized linear map `M` minimizing +``\|M \,\text{features} - \text{targets}\|^2 + \text{reg}\,\|M\|^2``. `features` is +``n_\text{in} \times T``, `targets` is ``n_\text{out} \times T``, and the returned +`M` has size ``n_\text{out} \times n_\text{in}``. +""" +function ridge_map(features::AbstractMatrix{<:Real}, targets::AbstractMatrix{<:Real}, reg::Real) + T = float(promote_type(eltype(features), eltype(targets))) + F = T.(features) + Y = T.(targets) + A = F * F' + T(reg) * I + B = F * Y' + return Matrix((A \ B)') +end + +# Per-pattern aperture lookup: a scalar applies to every pattern; a dictionary +# supplies a value per name. +aperture_for(aperture::Real, ::Symbol) = aperture +function aperture_for(aperture::AbstractDict{Symbol, <:Real}, name::Symbol) + haskey(aperture, name) || throw(ArgumentError("no aperture provided for :$name")) + return aperture[name] +end + +# ====================================================================== +# Loading patterns into a reservoir (Jaeger 2014, Section 3.2) +# ====================================================================== + +# Collect one pattern's reservoir states from a cleared carry, returning the +# post-washout states x(n), the lagged states x(n-1), and the aligned driver p(n). +function _drive_pattern( + rng::AbstractRNG, concept::Conceptor, signal::AbstractMatrix, + ps::NamedTuple, st::NamedTuple, washout::Int + ) + st_reset = resetcarry!(rng, concept.model, st.model; init_carry = nothing) + states, st_model = collectstates(concept.model, signal, ps.model, st_reset) + L = size(states, 2) + washout < L || throw(ArgumentError("washout=$washout ≥ signal length=$L")) + X = states[:, (washout + 1):L] # x(n) + Xlag = states[:, washout:(L - 1)] # x(n-1) + P = signal[:, (washout + 1):L] # driver aligned with x(n) + return X, Xlag, P, merge(st, (; model = st_model)) +end + +@doc raw""" + load!(rng, concept, named_signals, ps, st; + aperture=10.0, washout=500, reg_recurrent=1e-4, reg_readout=1e-2) + +Load named driving patterns into the reservoir wrapped by `concept` +(Jaeger 2014, Section 3.2). For each `name => signal` pair the reservoir is run +from a cleared carry, the first `washout` states are discarded, and the remaining +states are used to + + * derive and store a conceptor ``C = R (R + \alpha^{-2} I)^{-1}``, and + * accumulate data for a single shared input-internalizing recurrent matrix `W` + and readout `W_out`. + +`W` is fitted (ridge, `reg_recurrent`) so that +``W x(n-1) \approx W^* x(n-1) + W_\text{in} p(n)``, absorbing the input drive into +the recurrent weights; `W_out` is fitted (ridge, `reg_readout`) so that +``W_\text{out} x(n) \approx p(n)`` across all patterns. `aperture` is either a +scalar or a `name => α` dictionary. + +Returns `(ps, st)` with `reservoir_matrix` and `readout.weight` replaced and the +conceptor library populated. +""" +function load!( + rng::AbstractRNG, concept::Conceptor, named_signals, ps::NamedTuple, st::NamedTuple; + aperture = 10.0, washout::Int = 500, reg_recurrent::Real = 1.0e-4, reg_readout::Real = 1.0e-2 + ) + rps = reservoir_params(ps) + W_star = Float64.(rps.reservoir_matrix) + W_in = Float64.(rps.input_matrix) + + Xs = Matrix{Float64}[] + Xlags = Matrix{Float64}[] + Ps = Matrix{Float64}[] + st_acc = st + + for (name, signal) in named_signals + name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + sig = signal isa AbstractMatrix ? Float64.(signal) : reshape(Float64.(signal), 1, :) + X, Xlag, P, st_acc = _drive_pattern(rng, concept, sig, ps, st_acc, washout) + α = aperture_for(aperture, name) + store_conceptor!(st_acc, name, conceptor_from_states(X, α), α) + push!(Xs, X) + push!(Xlags, Xlag) + push!(Ps, P) + end + + Xall = reduce(hcat, Xs) + Xlagall = reduce(hcat, Xlags) + Pall = reduce(hcat, Ps) + + W_out = ridge_map(Xall, Pall, reg_readout) + recurrent_target = W_star * Xlagall .+ W_in * Pall + W = ridge_map(Xlagall, recurrent_target, reg_recurrent) + + T = eltype(rps.reservoir_matrix) + ps2 = set_reservoir_matrix(ps, T.(W)) + ps2 = set_readout_weight(ps2, T.(W_out)) + return ps2, st_acc +end + +# ====================================================================== +# Autonomous pattern generation (Jaeger 2014, Section 3.2) +# ====================================================================== + +@doc raw""" + generate(concept, ps, st; conceptor, steps, washout=200, + init_state=nothing, rng=Random.default_rng()) -> (Y, X) + +Run the loaded reservoir autonomously under `conceptor` and read out the observer +signal: + +```math +x(n) = C \tanh(W x(n-1) + b), \qquad y(n) = W_\text{out} x(n). +``` + +`conceptor` is either a stored conceptor `Symbol` or an explicit conceptor matrix +(e.g. the output of [`morph_conceptor`](@ref) or the Boolean operations). The first +`washout` steps let the autonomous orbit settle and are discarded. Returns `(Y, X)`: +the post-washout observer outputs (`out_dims × steps`) and reservoir states +(`res_dims × steps`). The rollout runs in `Float64` for numerical stability. +""" +function generate( + concept::Conceptor, ps::NamedTuple, st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix}, steps::Int, washout::Int = 200, + init_state::Union{Nothing, AbstractVector} = nothing, + rng::AbstractRNG = Random.default_rng() + ) + rps = reservoir_params(ps) + W = Float64.(rps.reservoir_matrix) + b = reservoir_bias(ps) + W_out = Float64.(readout_params(ps).weight) + C = resolve_conceptor(st, conceptor) + N = size(W, 1) + + x = init_state === nothing ? 0.5 .* randn(rng, N) : Float64.(init_state) + out_dims = size(W_out, 1) + Y = Matrix{Float64}(undef, out_dims, steps) + X = Matrix{Float64}(undef, N, steps) + + for n in 1:(washout + steps) + x = C * tanh.(W * x .+ b) + if n > washout + k = n - washout + @views X[:, k] = x + @views Y[:, k] = W_out * x + end + end + return Y, X +end + +# ====================================================================== +# Conceptor morphing (Jaeger 2014, Section 3.2) +# ====================================================================== + +@doc raw""" + morph_conceptor(st, weights) -> Matrix{Float64} + +Linearly combine stored conceptors into a morphed conceptor +``M = \sum_j \mu_j C_j`` (Jaeger 2014). `weights` is an iterable of `name => μ` +pairs, a `NamedTuple`, or a `Dict{Symbol,<:Real}`. Coefficients summing to one +interpolate between the named prototypes; coefficients outside `[0, 1]` extrapolate. +""" +function morph_conceptor(st::NamedTuple, weights) + M = nothing + for (name, μ) in pairs(weights) + C = get_conceptor(st, Symbol(name)) + C === nothing && throw(KeyError(Symbol(name))) + contribution = Float64(μ) .* C + M = M === nothing ? contribution : M .+ contribution + end + M === nothing && throw(ArgumentError("no weights provided")) + return M +end + +# ====================================================================== +# Conceptors as state filters for supervised readout training +# ====================================================================== + +# Normalize a target into a row-major matrix (1 × T for a vector target). +_as_matrix(t::AbstractMatrix) = t +_as_matrix(t::AbstractVector) = reshape(t, 1, :) + +""" + store_conceptors!(rng, concept, named_signals, ps, st; + aperture=1.0, init_carry=nothing) -> st + +Derive and store one conceptor per named signal without modifying the reservoir +weights. For each `name => signal` the reservoir is run from a cleared carry and +the conceptor `C(name) = R (R + α^{-2} I)^{-1}` of its state cloud is stored. Use +this to build a conceptor library for classification, where each class pattern +gets its own conceptor. `aperture` is a scalar or a `name => α` dictionary. +""" +function store_conceptors!( + rng::AbstractRNG, concept::Conceptor, named_signals, ps::NamedTuple, st::NamedTuple; + aperture = 1.0, init_carry = nothing + ) + st_acc = st + for (name, signal) in named_signals + name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + sig = signal isa AbstractMatrix ? signal : reshape(signal, 1, :) + st_reset = resetcarry!(rng, concept.model, st_acc.model; init_carry = init_carry) + states, st_model = collectstates(concept.model, sig, ps.model, st_reset) + st_acc = merge(st_acc, (; model = st_model)) + α = aperture_for(aperture, name) + store_conceptor!(st_acc, name, conceptor_from_states(states, α), α) + end + return st_acc +end + +""" + train!(rng, concept::Conceptor, named_signals, named_targets, ps, st, + train_method=StandardRidge(0.0); + washout=0, return_states=false, init_carry=nothing) -> (ps, st) + +Train a single readout on conceptor-filtered reservoir features. For each +`name => signal` the states are collected through the conceptor named `name` +(stored beforehand via [`store_conceptors!`](@ref)), the first `washout` columns +are dropped, and features and the matching `named_targets` are pooled across +patterns before the readout is fit with `train_method`. With `return_states=true` +the pooled feature matrix is also returned. +""" +function train!( + rng::AbstractRNG, concept::Conceptor, named_signals, named_targets, + ps::NamedTuple, st::NamedTuple, train_method = StandardRidge(0.0); + washout::Int = 0, return_states::Bool = false, init_carry = nothing, kwargs... + ) + targets = Dict{Symbol, AbstractMatrix}(Symbol(name) => _as_matrix(t) for (name, t) in named_targets) + + all_states = AbstractMatrix[] + all_targets = AbstractMatrix[] + for (name, signal) in named_signals + name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + haskey(targets, name) || throw(ArgumentError("no target data for pattern :$name")) + st_reset = resetcarry!(rng, concept.model, st.model; init_carry = init_carry) + st_tmp = merge(st, (; model = st_reset)) + states, _ = collectstates(concept, signal, ps, st_tmp; conceptor = name) + tgt = targets[name] + if washout > 0 + states = states[:, (washout + 1):end] + tgt = tgt[:, (washout + 1):end] + end + push!(all_states, states) + push!(all_targets, tgt) + end + + features = reduce(hcat, all_states) + pooled_targets = reduce(hcat, all_targets) + W_out = train(train_method, features, pooled_targets; kwargs...) + ps2, st2 = addreadout!(concept, W_out, ps, st) + return return_states ? ((ps2, st2), features) : (ps2, st2) +end diff --git a/test/test_conceptors.jl b/test/test_conceptors.jl new file mode 100644 index 00000000..f3fd426c --- /dev/null +++ b/test/test_conceptors.jl @@ -0,0 +1,131 @@ +using Test +using LinearAlgebra +using Random +using Statistics +using ReservoirComputing + +# A small reproducible reservoir state cloud for algebra tests. +function sample_correlation(rng, N, L) + X = randn(rng, N, L) + return correlation_matrix(X) +end + +@testset "Conceptors" begin + rng = Xoshiro(2024) + + @testset "conceptor matrix" begin + R = sample_correlation(rng, 20, 200) + C = conceptor_matrix(R, 10.0) + @test size(C) == (20, 20) + @test C ≈ C' # symmetric + ev = eigvals(Symmetric(C)) + @test all(>=(-1e-10), ev) # PSD + @test all(<=(1.0 + 1e-10), ev) # singular values ≤ 1 + @test conceptor_matrix(R, 1e8) ≈ I(20) atol = 1e-3 + @test norm(conceptor_matrix(R, 1e-6)) < 1e-6 + @test_throws ArgumentError conceptor_matrix(R, -1.0) + end + + @testset "aperture adaptation" begin + R = sample_correlation(rng, 15, 150) + C = conceptor_matrix(R, 5.0) + @test aperture_adapt(C, 1.0) ≈ C # γ = 1 is identity + # Prop 5: φ(φ(C, γ), β) = φ(C, γβ) + @test aperture_adapt(aperture_adapt(C, 2.0), 3.0) ≈ aperture_adapt(C, 6.0) + qs = [quota(aperture_adapt(C, g)) for g in (0.25, 1.0, 4.0, 16.0)] + @test issorted(qs) + P = aperture_adapt(C, Inf) + @test P * P ≈ P atol = 1e-8 # hardens to a projector + @test norm(aperture_adapt(C, 0.0)) < 1e-8 + @test conceptor_matrix(R, 7.0) ≈ aperture_adapt(conceptor_matrix(R, 1.0), 7.0) + end + + @testset "Boolean algebra" begin + C = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + B = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + A = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + + @test conceptor_not(conceptor_not(C)) ≈ C # double negation + @test conceptor_not(C) ≈ I(12) - C + # De Morgan + @test conceptor_or(C, B) ≈ conceptor_not(conceptor_and(conceptor_not(C), conceptor_not(B))) + @test conceptor_and(C, B) ≈ conceptor_not(conceptor_or(conceptor_not(C), conceptor_not(B))) + # commutativity + @test conceptor_and(C, B) ≈ conceptor_and(B, C) + @test conceptor_or(C, B) ≈ conceptor_or(B, C) + # AND shrinks, OR grows (quota ordering) + @test quota(conceptor_and(C, B)) <= quota(C) + 1e-8 + @test quota(C) <= quota(conceptor_or(C, B)) + 1e-8 + # associativity (holds for conceptors, unlike idempotence) + @test conceptor_and(conceptor_and(C, B), A) ≈ conceptor_and(C, conceptor_and(B, A)) atol = 1e-6 + @test conceptor_or(conceptor_or(C, B), A) ≈ conceptor_or(C, conceptor_or(B, A)) atol = 1e-6 + # neutral / absorbing elements + @test conceptor_and(C, Matrix(1.0I, 12, 12)) ≈ C atol = 1e-6 + @test conceptor_or(C, zeros(12, 12)) ≈ C atol = 1e-6 + # unexported infix sugar agrees + and_op, or_op, not_op = ReservoirComputing.:∧, ReservoirComputing.:∨, ReservoirComputing.:¬ + @test and_op(C, B) ≈ conceptor_and(C, B) + @test or_op(C, B) ≈ conceptor_or(C, B) + @test not_op(C) ≈ conceptor_not(C) + end + + @testset "wrapper / library" begin + rng2 = Xoshiro(7) + esn = ESN(1, 30, 1; use_bias = true) + concept = Conceptor(esn) + ps = initialparameters(rng2, concept) + st = initialstates(rng2, concept) + @test haskey(ps, :model) + @test isempty(st.conceptors) + C = conceptor_matrix(sample_correlation(rng2, 30, 100), 10.0) + store_conceptor!(st, :a, C, 10.0) + @test has_conceptor(st, :a) + @test get_conceptor(st, :a) ≈ C + st = set_active_conceptor(st, :a) + @test active_conceptor(st) ≈ C + @test_throws KeyError set_active_conceptor(st, :missing) + end + + @testset "load! / generate roundtrip" begin + rng3 = Xoshiro(99) + bias_init(r, dims...) = 0.2f0 .* randn(r, Float32, dims...) + input_init(r, dims...) = 1.5f0 .* randn(r, Float32, dims...) + res_init(r, dims...) = rand_sparse(r, dims...; radius = 1.5f0, sparsity = 0.1f0) + esn = ESN(1, 80, 1; use_bias = true, init_bias = bias_init, + init_input = input_init, init_reservoir = res_init) + concept = Conceptor(esn) + ps = initialparameters(rng3, concept) + st = initialstates(rng3, concept) + + n = 1:1000 + period = 8.8342522 + sig = reshape(Float32.(sin.(2pi .* n ./ period)), 1, :) + ps, st = load!(rng3, concept, [:sine => sig], ps, st; aperture = 1000.0, washout = 400) + @test has_conceptor(st, :sine) + + Y, X = generate(concept, ps, st; conceptor = :sine, steps = 200, washout = 300, rng = rng3) + @test size(Y) == (1, 200) + @test size(X) == (80, 200) + y = vec(Y) + best = minimum( + sqrt(mean(abs2, y .- [s * sin(2pi * k / period + ph) for k in eachindex(y)])) / std(y) + for ph in 0:0.02:2pi, s in (-1.0, 1.0) + ) + @test best < 0.05 + @test count(>(0.5), conceptor_singular_values(get_conceptor(st, :sine))) < 80 + end + + @testset "morph_conceptor" begin + rng4 = Xoshiro(5) + esn = ESN(1, 25, 1; use_bias = true) + concept = Conceptor(esn) + st = initialstates(rng4, concept) + Ca = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + Cb = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + store_conceptor!(st, :a, Ca, 10.0) + store_conceptor!(st, :b, Cb, 10.0) + @test morph_conceptor(st, (; a = 1.0, b = 0.0)) ≈ Ca + @test morph_conceptor(st, (; a = 0.3, b = 0.7)) ≈ 0.3 .* Ca .+ 0.7 .* Cb + @test_throws KeyError morph_conceptor(st, (; missing_name = 1.0)) + end +end From 6f03109f498487b4e1b9b227c82df859102be2fb Mon Sep 17 00:00:00 2001 From: MartinuzziFrancesco Date: Wed, 1 Jul 2026 16:18:41 +0200 Subject: [PATCH 2/2] feat: first pass of conceptors --- docs/src/assets/Project.toml | 41 ++ src/ReservoirComputing.jl | 2 +- src/conceptors.jl | 974 +++++++++++++++++++++++++---------- test/test_conceptors.jl | 200 ++++--- 4 files changed, 888 insertions(+), 329 deletions(-) create mode 100644 docs/src/assets/Project.toml diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml new file mode 100644 index 00000000..bdba183c --- /dev/null +++ b/docs/src/assets/Project.toml @@ -0,0 +1,41 @@ +[deps] +CellularAutomata = "878138dc-5b27-11ea-1a71-cb95d38d6b29" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" +DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LIBSVM = "b1bec4e5-fd48-53fe-b0cb-9723c09d164b" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +MLJLinearModels = "6ee0df7b-362f-4a72-a706-9e79364fb692" +OrdinaryDiffEqAdamsBashforthMoulton = "89bda076-bce5-4f1c-845f-551c83cdda9a" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReservoirComputing = "7c2d2b1e-3dd4-11ea-355a-8f6a8116e294" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" + +[compat] +CellularAutomata = "0.0.6" +ConcreteStructs = "0.2" +DataInterpolations = "6, 7, 8" +DelayDiffEq = "6" +Documenter = "1" +DocumenterCitations = "1" +DocumenterInterLinks = "1" +JLD2 = "0.6" +LIBSVM = "0.8" +LinearSolve = "3" +MLJLinearModels = "0.10" +OrdinaryDiffEqAdamsBashforthMoulton = "2" +OrdinaryDiffEqTsit5 = "1, 2" +Plots = "1" +ReservoirComputing = "0.12.0" +SciMLBase = "2.51, 3" +Static = "1" +StatsBase = "0.34.4" diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl index f38340a4..34962ccd 100644 --- a/src/ReservoirComputing.jl +++ b/src/ReservoirComputing.jl @@ -3,7 +3,7 @@ module ReservoirComputing using ArrayInterface: ArrayInterface using ConcreteStructs: @concrete using LinearAlgebra: eigvals, eigen, I, qr, Diagonal, diag, mul!, Symmetric, norm, - svd, svdvals, nullspace, pinv, tr + svd, svdvals, nullspace, pinv, tr, checksquare, issymmetric using LuxCore: AbstractLuxLayer, AbstractLuxContainerLayer, AbstractLuxWrapperLayer, setup, apply, replicate import LuxCore: initialparameters, initialstates, statelength, outputsize diff --git a/src/conceptors.jl b/src/conceptors.jl index 1c5edc17..75f8fbcf 100644 --- a/src/conceptors.jl +++ b/src/conceptors.jl @@ -1,79 +1,156 @@ -# Conceptors for recurrent neural networks, after Jaeger (2014), "Controlling -# Recurrent Neural Networks by Conceptors" (arXiv:1403.3369). -# -# A conceptor is a positive semidefinite matrix C = R (R + α^{-2} I)^{-1} derived -# from the correlation matrix R of a reservoir's driven states. It acts as a soft -# projector onto the linear subspace a driving pattern excites in the reservoir. -# Conceptors can be aperture-adapted, combined with a Boolean algebra, loaded into -# a reservoir to autonomously regenerate the loaded patterns, and morphed to -# interpolate or extrapolate between patterns. -# -# Conceptor matrices are always computed and stored in `Float64`: the defining -# operation is a matrix inversion whose conditioning degrades quickly for the -# near-singular correlation matrices produced by periodic drivers, so conceptor -# precision is decoupled from the (often `Float32`) reservoir precision. - -const ConceptorMatrix = Matrix{Float64} - -# ====================================================================== -# Core conceptor matrices -# ====================================================================== - @doc raw""" - correlation_matrix(states) -> Matrix{Float64} + correlation_matrix(states) -> Matrix State correlation matrix ``R = X X^\top / L`` for a reservoir state collection `states` of size ``N \times L`` (one state per column). This is the matrix whose regularized-identity map defines a conceptor. + +# Arguments + +- `states::AbstractMatrix{<:Real}`: Reservoir states arranged as + `(reservoir_dimension, sample_count)`. + +# Returns + +- A dense square matrix with the floating-point element type inferred from `states`. + +# Throws + +- `ArgumentError`: If `states` has no elements. + +# Example + +```julia +states = rand(Float32, 50, 1_000) +correlation = correlation_matrix(states) +``` """ function correlation_matrix(states::AbstractMatrix{<:Real}) - X = Float64.(states) - return (X * X') / size(X, 2) + isempty(states) && throw(ArgumentError("states must contain at least one state")) + element_type = float(eltype(states)) + state_matrix = element_type.(states) + return (state_matrix * state_matrix') / size(state_matrix, 2) end @doc raw""" - conceptor_matrix(R, aperture) -> Matrix{Float64} + conceptor_matrix(correlation, aperture) -> Matrix -Conceptor ``C = R (R + \alpha^{-2} I)^{-1}`` derived from a correlation matrix `R` -with aperture ``\alpha`` = `aperture` (Jaeger 2014, Eq. 7). The result is +Conceptor derived from a correlation matrix and the given `aperture` using +Jaeger (2014), Equation 7. The result is symmetric positive semidefinite with singular values in ``[0, 1)``. + +# Arguments + +- `correlation::AbstractMatrix{<:Real}`: Symmetric state correlation matrix. +- `aperture::Real`: Finite positive aperture. Larger values admit more state-space + directions; smaller values suppress more directions. + +# Returns + +- A dense conceptor matrix with the floating-point element type of `correlation`. + +# Throws + +- `ArgumentError`: If `aperture` is not finite and positive, or if `correlation` + is not symmetric. +- `DimensionMismatch`: If `correlation` is not square. + +# Example + +```julia +correlation = correlation_matrix(states) +conceptor = conceptor_matrix(correlation, 10) +``` """ -function conceptor_matrix(R::AbstractMatrix{<:Real}, aperture::Real) - aperture > 0 || throw(ArgumentError("aperture must be positive, got $aperture")) - Rd = Float64.(R) - inv_sq = 1.0 / aperture^2 - return Rd / (Rd + inv_sq * I) +function conceptor_matrix(correlation::AbstractMatrix{<:Real}, aperture::Real) + isfinite(aperture) && aperture > 0 || + throw(ArgumentError("aperture must be finite and positive, got $aperture")) + checksquare(correlation) + issymmetric(correlation) || throw(ArgumentError("correlation must be symmetric")) + element_type = float(eltype(correlation)) + correlation_matrix = element_type.(correlation) + typed_aperture = element_type(aperture) + inverse_aperture_squared = inv(typed_aperture * typed_aperture) + conceptor = if element_type <: Union{Float32, Float64} + decomposition = eigen(Symmetric(correlation_matrix)) + values = max.(decomposition.values, zero(element_type)) + decomposition.vectors * + Diagonal(values ./ (values .+ inverse_aperture_squared)) * + decomposition.vectors' + else + correlation_matrix / + Symmetric(correlation_matrix + inverse_aperture_squared * I) + end + return (conceptor + conceptor') / element_type(2) end """ - conceptor_from_states(states, aperture) -> Matrix{Float64} + conceptor_from_states(states, aperture) -> Matrix Convenience composition of [`correlation_matrix`](@ref) and [`conceptor_matrix`](@ref): the conceptor characterizing a reservoir state cloud `states` (size `N × L`) at the given `aperture`. + +# Arguments + +- `states::AbstractMatrix{<:Real}`: Reservoir states, one state per column. +- `aperture::Real`: Finite positive aperture. + +# Returns + +- The conceptor associated with the sample correlation of `states`. + +# Example + +```julia +conceptor = conceptor_from_states(states, 10) +``` """ function conceptor_from_states(states::AbstractMatrix{<:Real}, aperture::Real) return conceptor_matrix(correlation_matrix(states), aperture) end """ - conceptor_singular_values(C) -> Vector{Float64} + conceptor_singular_values(conceptor) -> Vector Singular values of a (symmetric) conceptor matrix `C`, in descending order. For a conceptor these coincide with its eigenvalues and lie in `[0, 1]`. + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- A vector of singular values in descending order. + +# Throws + +- `DimensionMismatch`: If `conceptor` is not square. """ -function conceptor_singular_values(C::AbstractMatrix{<:Real}) - return svdvals(Float64.(C)) +function conceptor_singular_values(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + return svdvals(float.(conceptor)) end """ - quota(C) -> Float64 + quota(conceptor) -> Real Mean singular value of a conceptor, `tr(C) / N`. Jaeger's "quota" measures the fraction of reservoir state space the conceptor leaves open (0 = point, 1 = all). + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- The mean singular value. Values near zero indicate a restrictive conceptor; + values near one indicate a permissive conceptor. """ -function quota(C::AbstractMatrix{<:Real}) - return tr(Float64.(C)) / size(C, 1) +function quota(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + return tr(conceptor) / size(conceptor, 1) end # ====================================================================== @@ -81,50 +158,116 @@ end # ====================================================================== @doc raw""" - adapt_singular_value(s, γ) -> Float64 + adapt_singular_value(singular_value, aperture_factor) -> Real -Aperture-adapted singular value ``s_\gamma`` for an input singular value -``s \in [0, 1]`` and adaptation factor ``\gamma \in [0, \infty]`` (Jaeger 2014, -Prop. 3, Eq. 19). Hard values `s = 0` and `s = 1` are fixed points; intermediate -values are pushed toward 0 as ``\gamma \to 0`` and toward 1 as ``\gamma \to \infty``. +Aperture-adapted value for an input `singular_value` in `[0, 1]` and a +nonnegative `aperture_factor` (Jaeger 2014, Proposition 3, Equation 19). +The boundary values zero and one are fixed points. + +# Arguments + +- `singular_value::Real`: Value to adapt, normally in `[0, 1]`. +- `aperture_factor::Real`: Nonnegative factor applied to the aperture. + +# Returns + +- The adapted singular value, using the floating-point type inferred from + `singular_value`. + +# Throws + +- `ArgumentError`: If `aperture_factor` is negative. """ -function adapt_singular_value(s::Float64, γ::Float64) - (s <= 0.0) && return 0.0 - (s >= 1.0) && return 1.0 - iszero(γ) && return 0.0 - isinf(γ) && return 1.0 - return s / (s + (1.0 - s) / γ^2) +function adapt_singular_value(singular_value::Real, aperture_factor::Real) + aperture_factor >= 0 || + throw(ArgumentError("aperture_factor must be nonnegative, got $aperture_factor")) + element_type = float(typeof(singular_value)) + typed_value = element_type(singular_value) + typed_factor = element_type(aperture_factor) + typed_value <= zero(element_type) && return zero(element_type) + typed_value >= one(element_type) && return one(element_type) + iszero(typed_factor) && return zero(element_type) + isinf(typed_factor) && return one(element_type) + return typed_value / + (typed_value + (one(element_type) - typed_value) / typed_factor^2) end @doc raw""" - aperture_adapt(C, γ) -> Matrix{Float64} + aperture_adapt(conceptor, aperture_factor) -> Matrix + +Aperture adaptation of a `conceptor` by a nonnegative `aperture_factor` +(Jaeger 2014, Definition 3). This is equivalent to multiplying its aperture by +the supplied factor. Zero and infinite factors are handled exactly. + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Symmetric square conceptor matrix. +- `aperture_factor::Real`: Nonnegative aperture multiplier. + +# Returns -Aperture adaptation ``\varphi(C, \gamma)`` of a conceptor `C` by factor -``\gamma \in [0, \infty]`` (Jaeger 2014, Def. 3). Equivalent to multiplying the -aperture of `C` by ``\gamma``; in particular ``C(R, \alpha) = \varphi(C(R, 1), \alpha)``. -Applied via the eigendecomposition so the limiting cases ``\gamma = 0`` (hard -zero) and ``\gamma = \infty`` (hardening to a projector) are exact. +- A conceptor with the same principal directions and adapted singular values. + +# Throws + +- `ArgumentError`: If the factor is negative or `conceptor` is not symmetric. +- `DimensionMismatch`: If `conceptor` is not square. + +# Example + +```julia +wider_conceptor = aperture_adapt(conceptor, 2) +projector = aperture_adapt(conceptor, Inf) +``` """ -function aperture_adapt(C::AbstractMatrix{<:Real}, γ::Real) - γ >= 0 || throw(ArgumentError("aperture factor γ must be ≥ 0, got $γ")) - F = eigen(Symmetric(Float64.(C))) - s_adapted = adapt_singular_value.(clamp.(F.values, 0.0, 1.0), Float64(γ)) - return F.vectors * Diagonal(s_adapted) * F.vectors' +function aperture_adapt(conceptor::AbstractMatrix{<:Real}, aperture_factor::Real) + aperture_factor >= 0 || + throw(ArgumentError("aperture_factor must be nonnegative, got $aperture_factor")) + checksquare(conceptor) + issymmetric(conceptor) || throw(ArgumentError("conceptor must be symmetric")) + element_type = float(eltype(conceptor)) + conceptor_matrix = element_type.(conceptor) + typed_factor = element_type(aperture_factor) + if isfinite(typed_factor) && !iszero(typed_factor) + inverse_factor_squared = inv(typed_factor * typed_factor) + adapted = conceptor_matrix / + Symmetric(conceptor_matrix + inverse_factor_squared * (I - conceptor_matrix)) + return (adapted + adapted') / element_type(2) + end + decomposition = eigen(Symmetric(conceptor_matrix)) + adapted_values = adapt_singular_value.( + clamp.(decomposition.values, zero(element_type), one(element_type)), typed_factor + ) + return decomposition.vectors * Diagonal(adapted_values) * decomposition.vectors' end """ - reaperture(C, from_aperture, to_aperture) -> Matrix{Float64} + reaperture(conceptor, from_aperture, to_aperture) -> Matrix Re-express a conceptor `C` that currently has aperture `from_aperture` so that it -has aperture `to_aperture`, using `φ(C, to/from)` (Jaeger 2014, Prop. 5). +has aperture `to_aperture` (Jaeger 2014, Proposition 5). + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Conceptor formed at `from_aperture`. +- `from_aperture::Real`: Finite positive current aperture. +- `to_aperture::Real`: Nonnegative target aperture; `Inf` is allowed. + +# Returns + +- The aperture-adapted conceptor. """ -function reaperture(C::AbstractMatrix{<:Real}, from_aperture::Real, to_aperture::Real) - from_aperture > 0 || throw(ArgumentError("from_aperture must be positive")) - return aperture_adapt(C, to_aperture / from_aperture) +function reaperture( + conceptor::AbstractMatrix{<:Real}, from_aperture::Real, to_aperture::Real + ) + isfinite(from_aperture) && from_aperture > 0 || + throw(ArgumentError("from_aperture must be finite and positive")) + to_aperture >= 0 || throw(ArgumentError("to_aperture must be nonnegative")) + return aperture_adapt(conceptor, to_aperture / from_aperture) end @doc raw""" - attenuation(C, W, b; steps, washout, init_state, rng) -> Float64 + attenuation(conceptor, recurrent_weights, bias; steps, washout, init_state, rng) -> Real Attenuation ``a_C = E[\|z(n) - x(n)\|^2] / E[\|z(n)\|^2]`` of a loaded reservoir (recurrent weights `W`, bias `b`) run autonomously under conceptor `C` @@ -132,48 +275,106 @@ Attenuation ``a_C = E[\|z(n) - x(n)\|^2] / E[\|z(n)\|^2]`` of a loaded reservoir update and ``x(n) = C z(n)`` is the conceptor-constrained state. The attenuation is the fraction of reservoir signal energy suppressed by `C`; as a function of aperture it passes through a minimum at the best-reconstructing aperture. + +# Arguments + +- `conceptor`: Square conceptor matching the recurrent reservoir dimension. +- `recurrent_weights`: Square autonomous recurrent weight matrix. +- `bias`: Reservoir bias vector. + +# Keywords + +- `steps::Int = 500`: Number of samples included in the energy ratio. +- `washout::Int = 200`: Initial autonomous steps excluded from the ratio. +- `init_state = nothing`: Optional initial reservoir state. +- `rng = Random.default_rng()`: Random number generator used when `init_state` is + omitted. + +# Returns + +- The fraction of unconstrained reservoir energy suppressed by the conceptor. """ function attenuation( - C::AbstractMatrix{<:Real}, W::AbstractMatrix{<:Real}, b::AbstractVector{<:Real}; - steps::Int = 500, washout::Int = 200, + conceptor::AbstractMatrix{<:Real}, + recurrent_weights::AbstractMatrix{<:Real}, + bias::AbstractVector{<:Real}; + steps::Int = 500, + washout::Int = 200, init_state::Union{Nothing, AbstractVector} = nothing, - rng::AbstractRNG = Random.default_rng() + rng::AbstractRNG = Random.default_rng(), ) - Cd = Float64.(C) - Wd = Float64.(W) - bd = Float64.(b) - N = size(Wd, 1) - x = init_state === nothing ? 0.5 .* randn(rng, N) : Float64.(init_state) - - num = 0.0 - den = 0.0 - for n in 1:(washout + steps) - z = tanh.(Wd * x .+ bd) - x = Cd * z - if n > washout - num += sum(abs2, z .- x) - den += sum(abs2, z) + checksquare(conceptor) + checksquare(recurrent_weights) + size(conceptor) == size(recurrent_weights) || + throw(DimensionMismatch("conceptor and recurrent_weights must have equal dimensions")) + size(recurrent_weights, 1) == length(bias) || + throw(DimensionMismatch("bias length must equal the recurrent matrix dimension")) + steps > 0 || throw(ArgumentError("steps must be positive")) + washout >= 0 || throw(ArgumentError("washout must be nonnegative")) + element_type = float( + promote_type(eltype(conceptor), eltype(recurrent_weights), eltype(bias)) + ) + conceptor_matrix = element_type.(conceptor) + recurrent_matrix = element_type.(recurrent_weights) + bias_vector = element_type.(bias) + reservoir_dimension = size(recurrent_matrix, 1) + state = init_state === nothing ? + element_type(0.5) .* randn(rng, element_type, reservoir_dimension) : + element_type.(init_state) + length(state) == reservoir_dimension || + throw(DimensionMismatch("init_state must have length $reservoir_dimension")) + + suppressed_energy = zero(element_type) + total_energy = zero(element_type) + for step in 1:(washout + steps) + unconstrained_state = tanh.(recurrent_matrix * state .+ bias_vector) + state = conceptor_matrix * unconstrained_state + if step > washout + suppressed_energy += sum(abs2, unconstrained_state .- state) + total_energy += sum(abs2, unconstrained_state) end end - return num / den + iszero(total_energy) && + throw(ArgumentError("attenuation is undefined for a zero-energy rollout")) + return suppressed_energy / total_energy end """ - optimal_aperture(R, W, b, apertures; kwargs...) -> (best_aperture, attenuations) + optimal_aperture(correlation, recurrent_weights, bias, apertures; kwargs...) Select the aperture minimizing the [`attenuation`](@ref) criterion over a grid -`apertures`. For each candidate `α` the conceptor `C(R, α)` is formed and the +`apertures`. For each candidate aperture, its conceptor is formed and the loaded reservoir (`W`, `b`) is run autonomously to measure its attenuation. Returns the minimizing aperture together with the full vector of attenuations (aligned with `apertures`) so the characteristic trough can be inspected. `kwargs` are forwarded to [`attenuation`](@ref). + +# Arguments + +- `correlation`: Symmetric state correlation matrix. +- `recurrent_weights`: Loaded recurrent weight matrix. +- `bias`: Reservoir bias vector. +- `apertures`: Nonempty vector of finite positive candidate apertures. + +# Returns + +- `(best_aperture, attenuations)`, where `attenuations` follows the order of + `apertures`. """ function optimal_aperture( - R::AbstractMatrix{<:Real}, W::AbstractMatrix{<:Real}, b::AbstractVector{<:Real}, - apertures::AbstractVector{<:Real}; kwargs... + correlation::AbstractMatrix{<:Real}, + recurrent_weights::AbstractMatrix{<:Real}, + bias::AbstractVector{<:Real}, + apertures::AbstractVector{<:Real}; + kwargs..., ) - atts = [attenuation(conceptor_matrix(R, α), W, b; kwargs...) for α in apertures] - return apertures[argmin(atts)], atts + isempty(apertures) && throw(ArgumentError("apertures must not be empty")) + attenuations = [ + attenuation( + conceptor_matrix(correlation, aperture), recurrent_weights, bias; kwargs... + ) for aperture in apertures + ] + return apertures[argmin(attenuations)], attenuations end # ====================================================================== @@ -181,71 +382,110 @@ end # ====================================================================== @doc raw""" - conceptor_not(C) -> Matrix{Float64} + conceptor_not(conceptor) -> Matrix -Logical negation ``\neg C = I - C`` of a conceptor (Jaeger 2014, Eq. 28). +Logical negation of a conceptor (Jaeger 2014, Equation 28). Exchanges the roles of the directions the conceptor admits and suppresses. -""" -function conceptor_not(C::AbstractMatrix{<:Real}) - Cd = Float64.(C) - return Matrix(I, size(Cd)) - Cd -end -# Orthonormal basis of the range of a symmetric PSD matrix, using a tolerance on -# the singular values to decide the numerical rank. -function _range_basis(C::Matrix{Float64}; tol::Float64 = 1.0e-12) - F = svd(C) - rank = count(>(tol), F.S) - return F.U[:, 1:rank] +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- The complementary conceptor. Its singular values are one minus those of + `conceptor`. +""" +function conceptor_not(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + element_type = float(eltype(conceptor)) + conceptor_matrix = element_type.(conceptor) + return Matrix{element_type}(I, size(conceptor_matrix)) - conceptor_matrix end @doc raw""" - conceptor_and(C, B) -> Matrix{Float64} + conceptor_and(first_conceptor, second_conceptor) -> Matrix -Logical conjunction ``C \wedge B`` of two conceptors (Jaeger 2014, Eq. 32). For +Logical conjunction of two conceptors (Jaeger 2014, Equation 32). For full-rank conceptors this equals ``(C^{-1} + B^{-1} - I)^{-1}``; the implementation uses the range-intersection projector form so that singular conceptors are handled robustly. The result admits exactly the reservoir directions admitted by *both* `C` and `B`. + +# Arguments + +- `first_conceptor::AbstractMatrix{<:Real}`: First symmetric conceptor. +- `second_conceptor::AbstractMatrix{<:Real}`: Second symmetric conceptor of the + same size. + +# Returns + +- The conjunction of the two conceptors, including when either input is + rank-deficient. + +# Throws + +- `DimensionMismatch`: If the conceptors are not square or have different sizes. +- `ArgumentError`: If either conceptor is not symmetric. """ -function conceptor_and(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) - Cd = Float64.(C) - Bd = Float64.(B) - size(Cd) == size(Bd) || throw(DimensionMismatch("conceptors must have equal size")) +function conceptor_and( + first_conceptor::AbstractMatrix{<:Real}, + second_conceptor::AbstractMatrix{<:Real} + ) + checksquare(first_conceptor) + checksquare(second_conceptor) + element_type = float(promote_type(eltype(first_conceptor), eltype(second_conceptor))) + first_matrix = element_type.(first_conceptor) + second_matrix = element_type.(second_conceptor) + size(first_matrix) == size(second_matrix) || + throw(DimensionMismatch("conceptors must have equal size")) + issymmetric(first_matrix) || throw(ArgumentError("first_conceptor must be symmetric")) + issymmetric(second_matrix) || throw(ArgumentError("second_conceptor must be symmetric")) # Basis of R(C) ∩ R(B) = N(P_{N(C)} + P_{N(B)}), via the null-space projectors. - UC0 = nullspace(Cd) - UB0 = nullspace(Bd) - M = UC0 * UC0' + UB0 * UB0' - Fm = svd(M) - rank_M = count(>(1.0e-12), Fm.S) - Bgk = Fm.U[:, (rank_M + 1):end] # basis of the range intersection - - if size(Bgk, 2) == 0 - return zeros(size(Cd)) # disjoint supports → all-zero conceptor + first_nullspace = nullspace(first_matrix) + second_nullspace = nullspace(second_matrix) + nullspace_projector = + first_nullspace * first_nullspace' + second_nullspace * second_nullspace' + intersection_basis = nullspace(nullspace_projector) + + if size(intersection_basis, 2) == 0 + return zeros(element_type, size(first_matrix)) end - core = Bgk' * (pinv(Cd) + pinv(Bd) - I) * Bgk - return Bgk * (core \ Matrix(I, size(core))) * Bgk' + core = intersection_basis' * (pinv(first_matrix) + pinv(second_matrix) - I) * + intersection_basis + result = intersection_basis * + (core \ Matrix{element_type}(I, size(core))) * + intersection_basis' + return (result + result') / element_type(2) end @doc raw""" - conceptor_or(C, B) -> Matrix{Float64} + conceptor_or(first_conceptor, second_conceptor) -> Matrix -Logical disjunction ``C \vee B`` of two conceptors (Jaeger 2014, Eq. 31), defined -via De Morgan's law ``C \vee B = \neg(\neg C \wedge \neg B)``. The result admits +Logical disjunction of two conceptors (Jaeger 2014, Equation 31), defined +via De Morgan's law. The result admits every reservoir direction admitted by `C` or by `B`. + +# Arguments + +- `first_conceptor::AbstractMatrix{<:Real}`: First symmetric conceptor. +- `second_conceptor::AbstractMatrix{<:Real}`: Second symmetric conceptor of the + same size. + +# Returns + +- The disjunction of the two conceptors. """ -function conceptor_or(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) - return conceptor_not(conceptor_and(conceptor_not(C), conceptor_not(B))) +function conceptor_or( + first_conceptor::AbstractMatrix{<:Real}, + second_conceptor::AbstractMatrix{<:Real} + ) + return conceptor_not( + conceptor_and(conceptor_not(first_conceptor), conceptor_not(second_conceptor)) + ) end -# Infix sugar mirroring the paper's ¬ / ∧ / ∨ notation. `∧` and `∨` are Julia infix -# operators; `¬` is a prefix-style function (`¬(C)`). Not exported, to avoid -# polluting the package namespace; use the named functions or `ReservoirComputing.:∧`. -const ¬ = conceptor_not -∧(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) = conceptor_and(C, B) -∨(C::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real}) = conceptor_or(C, B) - # ====================================================================== # The Conceptor reservoir-computer wrapper and its parameter/state plumbing # ====================================================================== @@ -260,9 +500,27 @@ algebra, and used to constrain autonomous generation (Jaeger 2014). The wrapped `model` provides the reservoir update; the `Conceptor` adds a conceptor library to the model state. Parameters and states are nested under the `model` field, leaving room for conceptor bookkeeping at the top level. + +# Arguments + +- `model`: Reservoir computer to wrap. Pattern loading expects an ESN-compatible + parameter layout with reservoir, input, bias, and readout weights. + +# Returns + +- A conceptor-enabled reservoir computer used with [`initialparameters`](@ref), + [`initialstates`](@ref), [`load!`](@ref), and [`generate`](@ref). + +# Example + +```julia +model = Conceptor(ESN(1, 100, 1; use_bias = true)) +parameters = initialparameters(rng, model) +states = initialstates(rng, model) +``` """ @concrete struct Conceptor <: AbstractReservoirComputer{(:model,)} - model + model::Any end function initialparameters(rng::AbstractRNG, concept::Conceptor) @@ -272,8 +530,8 @@ end function initialstates(rng::AbstractRNG, concept::Conceptor) return (; model = initialstates(rng, concept.model), - conceptors = Dict{Symbol, ConceptorMatrix}(), - apertures = Dict{Symbol, Float64}(), + conceptors = Dict{Symbol, Matrix}(), + apertures = Dict{Symbol, Real}(), active_conceptor = nothing, ) end @@ -282,25 +540,48 @@ end has_conceptor(st, name) -> Bool Whether a conceptor called `name` has been stored in the state `st`. + +Returns `true` when `name` is present and `false` otherwise. """ has_conceptor(st::NamedTuple, name::Symbol) = haskey(st.conceptors, name) """ - get_conceptor(st, name) -> Matrix{Float64} or nothing + get_conceptor(st, name) -> Matrix or nothing The stored conceptor matrix called `name`, or `nothing` if it is absent. + +Use [`has_conceptor`](@ref) when absence should be checked separately. """ get_conceptor(st::NamedTuple, name::Symbol) = get(st.conceptors, name, nothing) """ - store_conceptor!(st, name, C, aperture) -> st + store_conceptor!(st, name, conceptor, aperture) -> st Record conceptor matrix `C` (formed at `aperture`) under `name` in the state's conceptor library. Mutates the library dictionaries in place and returns `st`. + +# Arguments + +- `st::NamedTuple`: State returned by `initialstates` for a [`Conceptor`](@ref). +- `name::Symbol`: Name used to retrieve the conceptor later. +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. +- `aperture::Real`: Finite positive aperture associated with the matrix. + +# Returns + +- The same state object, with its conceptor and aperture dictionaries updated. """ -function store_conceptor!(st::NamedTuple, name::Symbol, C::AbstractMatrix{<:Real}, aperture::Real) - st.conceptors[name] = Float64.(C) - st.apertures[name] = Float64(aperture) +function store_conceptor!( + st::NamedTuple, + name::Symbol, + conceptor::AbstractMatrix{<:Real}, + aperture::Real, + ) + checksquare(conceptor) + isfinite(aperture) && aperture > 0 || + throw(ArgumentError("aperture must be finite and positive, got $aperture")) + st.conceptors[name] = Matrix(conceptor) + st.apertures[name] = convert(float(eltype(conceptor)), aperture) return st end @@ -309,6 +590,9 @@ end Mark the conceptor `name` (or `nothing`) as the active one, returning an updated state. The active conceptor is the default constraint for generation. + +Passing `nothing` clears the active selection. A `KeyError` is thrown when a +nonexistent name is selected. """ function set_active_conceptor(st::NamedTuple, name::Union{Symbol, Nothing}) name !== nothing && !has_conceptor(st, name) && throw(KeyError(name)) @@ -316,22 +600,25 @@ function set_active_conceptor(st::NamedTuple, name::Union{Symbol, Nothing}) end """ - active_conceptor(st) -> Matrix{Float64} + active_conceptor(st) -> Matrix The currently active conceptor matrix. Throws if none has been set. + +Use [`set_active_conceptor`](@ref) to select or clear the active conceptor. """ function active_conceptor(st::NamedTuple) name = st.active_conceptor - name === nothing && throw(ArgumentError("no active conceptor; call set_active_conceptor")) + name === nothing && + throw(ArgumentError("no active conceptor; call set_active_conceptor")) return st.conceptors[name] end # Resolve a conceptor argument that is either an explicit matrix or a stored name. -resolve_conceptor(st::NamedTuple, conceptor::AbstractMatrix{<:Real}) = Float64.(conceptor) +resolve_conceptor(::NamedTuple, conceptor::AbstractMatrix{<:Real}) = conceptor function resolve_conceptor(st::NamedTuple, conceptor::Symbol) - C = get_conceptor(st, conceptor) - C === nothing && throw(KeyError(conceptor)) - return C + stored_conceptor = get_conceptor(st, conceptor) + stored_conceptor === nothing && throw(KeyError(conceptor)) + return stored_conceptor end """ @@ -343,8 +630,11 @@ states are filtered through it, `C x(n)`, realizing the conceptor as a state modifier. """ function collectstates( - concept::Conceptor, signal::AbstractMatrix, ps::NamedTuple, st::NamedTuple; - conceptor::Union{Symbol, AbstractMatrix, Nothing} = nothing + concept::Conceptor, + signal::AbstractMatrix, + ps::NamedTuple, + st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix, Nothing} = nothing, ) states, st_model = collectstates(concept.model, signal, ps.model, st.model) new_st = merge(st, (; model = st_model)) @@ -353,35 +643,35 @@ function collectstates( end """ - addreadout!(concept::Conceptor, W, ps, st) -> (ps, st) + addreadout!(concept::Conceptor, weights, ps, st) -> (ps, st) Install trained readout weights `W` into the wrapped model's readout layer (`ps.model.readout.weight`), keeping the generation and classification paths on a single readout convention. """ -function addreadout!(concept::Conceptor, W::AbstractMatrix, ps::NamedTuple, st::NamedTuple) - return set_readout_weight(ps, W), st +function addreadout!( + concept::Conceptor, weights::AbstractMatrix, ps::NamedTuple, st::NamedTuple + ) + return set_readout_weight(ps, weights), st end -# ====================================================================== -# Parameter-tree access for the wrapped model (ps.model.{reservoir,readout}) -# ====================================================================== - reservoir_params(ps::NamedTuple) = ps.model.reservoir readout_params(ps::NamedTuple) = ps.model.readout function reservoir_bias(ps::NamedTuple) rps = reservoir_params(ps) - return haskey(rps, :bias) ? Float64.(vec(rps.bias)) : zeros(Float64, size(rps.reservoir_matrix, 1)) + element_type = eltype(rps.reservoir_matrix) + return haskey(rps, :bias) ? + vec(rps.bias) : zeros(element_type, size(rps.reservoir_matrix, 1)) end -function set_reservoir_matrix(ps::NamedTuple, W::AbstractMatrix) - reservoir = merge(ps.model.reservoir, (; reservoir_matrix = W)) +function set_reservoir_matrix(ps::NamedTuple, weights::AbstractMatrix) + reservoir = merge(ps.model.reservoir, (; reservoir_matrix = weights)) return merge(ps, (; model = merge(ps.model, (; reservoir = reservoir)))) end -function set_readout_weight(ps::NamedTuple, W::AbstractMatrix) - readout = merge(ps.model.readout, (; weight = W)) +function set_readout_weight(ps::NamedTuple, weights::AbstractMatrix) + readout = merge(ps.model.readout, (; weight = weights)) return merge(ps, (; model = merge(ps.model, (; readout = readout)))) end @@ -392,14 +682,26 @@ Tikhonov-regularized linear map `M` minimizing ``\|M \,\text{features} - \text{targets}\|^2 + \text{reg}\,\|M\|^2``. `features` is ``n_\text{in} \times T``, `targets` is ``n_\text{out} \times T``, and the returned `M` has size ``n_\text{out} \times n_\text{in}``. + +# Arguments + +- `features`: Feature matrix with observations in columns. +- `targets`: Target matrix with observations in columns. +- `reg::Real`: Nonnegative ridge penalty. + +# Returns + +- A matrix mapping feature columns to target columns. """ -function ridge_map(features::AbstractMatrix{<:Real}, targets::AbstractMatrix{<:Real}, reg::Real) - T = float(promote_type(eltype(features), eltype(targets))) - F = T.(features) - Y = T.(targets) - A = F * F' + T(reg) * I - B = F * Y' - return Matrix((A \ B)') +function ridge_map( + features::AbstractMatrix{<:Real}, + targets::AbstractMatrix{<:Real}, + reg::Real, + ) + element_type = float(promote_type(eltype(features), eltype(targets))) + feature_matrix = element_type.(features) + target_matrix = element_type.(targets) + return train(StandardRidge(element_type, reg), feature_matrix, target_matrix) end # Per-pattern aperture lookup: a scalar applies to every pattern; a dictionary @@ -410,24 +712,25 @@ function aperture_for(aperture::AbstractDict{Symbol, <:Real}, name::Symbol) return aperture[name] end -# ====================================================================== -# Loading patterns into a reservoir (Jaeger 2014, Section 3.2) -# ====================================================================== +# loading patterns -# Collect one pattern's reservoir states from a cleared carry, returning the -# post-washout states x(n), the lagged states x(n-1), and the aligned driver p(n). function _drive_pattern( - rng::AbstractRNG, concept::Conceptor, signal::AbstractMatrix, - ps::NamedTuple, st::NamedTuple, washout::Int + rng::AbstractRNG, + concept::Conceptor, + signal::AbstractMatrix, + ps::NamedTuple, + st::NamedTuple, + washout::Int, ) st_reset = resetcarry!(rng, concept.model, st.model; init_carry = nothing) states, st_model = collectstates(concept.model, signal, ps.model, st_reset) - L = size(states, 2) - washout < L || throw(ArgumentError("washout=$washout ≥ signal length=$L")) - X = states[:, (washout + 1):L] # x(n) - Xlag = states[:, washout:(L - 1)] # x(n-1) - P = signal[:, (washout + 1):L] # driver aligned with x(n) - return X, Xlag, P, merge(st, (; model = st_model)) + signal_length = size(states, 2) + washout < signal_length || + throw(ArgumentError("washout=$washout exceeds signal length=$signal_length")) + current_states = states[:, (washout + 1):signal_length] + lagged_states = states[:, washout:(signal_length - 1)] + aligned_signal = signal[:, (washout + 1):signal_length] + return current_states, lagged_states, aligned_signal, merge(st, (; model = st_model)) end @doc raw""" @@ -439,7 +742,7 @@ Load named driving patterns into the reservoir wrapped by `concept` from a cleared carry, the first `washout` states are discarded, and the remaining states are used to - * derive and store a conceptor ``C = R (R + \alpha^{-2} I)^{-1}``, and + * derive and store a conceptor from the state correlation matrix, and * accumulate data for a single shared input-internalizing recurrent matrix `W` and readout `W_out`. @@ -447,47 +750,89 @@ states are used to ``W x(n-1) \approx W^* x(n-1) + W_\text{in} p(n)``, absorbing the input drive into the recurrent weights; `W_out` is fitted (ridge, `reg_readout`) so that ``W_\text{out} x(n) \approx p(n)`` across all patterns. `aperture` is either a -scalar or a `name => α` dictionary. +scalar or a dictionary mapping each name to an aperture. Returns `(ps, st)` with `reservoir_matrix` and `readout.weight` replaced and the conceptor library populated. + +# Arguments + +- `rng::AbstractRNG`: Random number generator used when reservoir carries reset. +- `concept::Conceptor`: Conceptor-wrapped reservoir. +- `named_signals`: Iterable of `Symbol => signal` pairs. A vector signal is treated + as one-dimensional; a matrix has one input channel per row and time per column. +- `ps::NamedTuple`: Current parameters. +- `st::NamedTuple`: Current states. + +# Keywords + +- `aperture = 10.0`: One aperture for every pattern, or a dictionary keyed by + pattern name. +- `washout::Int = 500`: Initial samples excluded from fitting. +- `reg_recurrent::Real = 1.0e-4`: Ridge penalty for recurrent weights. +- `reg_readout::Real = 1.0e-2`: Ridge penalty for readout weights. + +# Returns + +- `(parameters, states)`: Updated parameters and a state containing one conceptor + per named signal. + +# Throws + +- `ArgumentError`: If a name is not a `Symbol`, an aperture is missing, or washout + is not shorter than a signal. """ function load!( - rng::AbstractRNG, concept::Conceptor, named_signals, ps::NamedTuple, st::NamedTuple; - aperture = 10.0, washout::Int = 500, reg_recurrent::Real = 1.0e-4, reg_readout::Real = 1.0e-2 + rng::AbstractRNG, + concept::Conceptor, + named_signals, + ps::NamedTuple, + st::NamedTuple; + aperture = 10.0, + washout::Int = 500, + reg_recurrent::Real = 1.0e-4, + reg_readout::Real = 1.0e-2, ) rps = reservoir_params(ps) - W_star = Float64.(rps.reservoir_matrix) - W_in = Float64.(rps.input_matrix) + element_type = float(eltype(rps.reservoir_matrix)) + original_recurrent_weights = element_type.(rps.reservoir_matrix) + input_weights = element_type.(rps.input_matrix) - Xs = Matrix{Float64}[] - Xlags = Matrix{Float64}[] - Ps = Matrix{Float64}[] + state_sets = Matrix{element_type}[] + lagged_state_sets = Matrix{element_type}[] + signal_sets = Matrix{element_type}[] st_acc = st for (name, signal) in named_signals - name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) - sig = signal isa AbstractMatrix ? Float64.(signal) : reshape(Float64.(signal), 1, :) - X, Xlag, P, st_acc = _drive_pattern(rng, concept, sig, ps, st_acc, washout) - α = aperture_for(aperture, name) - store_conceptor!(st_acc, name, conceptor_from_states(X, α), α) - push!(Xs, X) - push!(Xlags, Xlag) - push!(Ps, P) + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + signal_matrix = signal isa AbstractMatrix ? + element_type.(signal) : reshape(element_type.(signal), 1, :) + current_states, lagged_states, aligned_signal, st_acc = + _drive_pattern(rng, concept, signal_matrix, ps, st_acc, washout) + pattern_aperture = aperture_for(aperture, name) + store_conceptor!( + st_acc, name, conceptor_from_states(current_states, pattern_aperture), + pattern_aperture + ) + push!(state_sets, current_states) + push!(lagged_state_sets, lagged_states) + push!(signal_sets, aligned_signal) end - Xall = reduce(hcat, Xs) - Xlagall = reduce(hcat, Xlags) - Pall = reduce(hcat, Ps) + all_states = reduce(hcat, state_sets) + all_lagged_states = reduce(hcat, lagged_state_sets) + all_signals = reduce(hcat, signal_sets) - W_out = ridge_map(Xall, Pall, reg_readout) - recurrent_target = W_star * Xlagall .+ W_in * Pall - W = ridge_map(Xlagall, recurrent_target, reg_recurrent) + readout_weights = ridge_map(all_states, all_signals, reg_readout) + recurrent_target = + original_recurrent_weights * all_lagged_states .+ input_weights * all_signals + recurrent_weights = ridge_map(all_lagged_states, recurrent_target, reg_recurrent) - T = eltype(rps.reservoir_matrix) - ps2 = set_reservoir_matrix(ps, T.(W)) - ps2 = set_readout_weight(ps2, T.(W_out)) - return ps2, st_acc + updated_parameters = set_reservoir_matrix(ps, element_type.(recurrent_weights)) + updated_parameters = + set_readout_weight(updated_parameters, element_type.(readout_weights)) + return updated_parameters, st_acc end # ====================================================================== @@ -509,59 +854,126 @@ x(n) = C \tanh(W x(n-1) + b), \qquad y(n) = W_\text{out} x(n). (e.g. the output of [`morph_conceptor`](@ref) or the Boolean operations). The first `washout` steps let the autonomous orbit settle and are discarded. Returns `(Y, X)`: the post-washout observer outputs (`out_dims × steps`) and reservoir states -(`res_dims × steps`). The rollout runs in `Float64` for numerical stability. +(`res_dims × steps`). The rollout uses the element type of the reservoir parameters. + +# Arguments + +- `concept::Conceptor`: Loaded conceptor model. +- `ps::NamedTuple`: Parameters returned by [`load!`](@ref). +- `st::NamedTuple`: State containing the requested conceptor. + +# Keywords + +- `conceptor`: Stored conceptor name or an explicit square conceptor matrix. +- `steps::Int`: Number of returned time steps; must be positive. +- `washout::Int = 200`: Autonomous steps discarded before recording. +- `init_state = nothing`: Optional reservoir state vector. A random state is used + when omitted. +- `rng = Random.default_rng()`: Random number generator for the initial state. + +# Returns + +- `(outputs, states)`: Observer outputs and reservoir states, with time along + columns. + +# Throws + +- `KeyError`: If a requested conceptor name is absent. +- `DimensionMismatch`: If the conceptor or initial state does not match the + reservoir dimension. +- `ArgumentError`: If `steps` is not positive or `washout` is negative. + +# Example + +```julia +outputs, states = generate( + concept, parameters, model_state; + conceptor = :sine, steps = 500, washout = 100, rng +) +``` """ function generate( - concept::Conceptor, ps::NamedTuple, st::NamedTuple; - conceptor::Union{Symbol, AbstractMatrix}, steps::Int, washout::Int = 200, + concept::Conceptor, + ps::NamedTuple, + st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix}, + steps::Int, + washout::Int = 200, init_state::Union{Nothing, AbstractVector} = nothing, - rng::AbstractRNG = Random.default_rng() + rng::AbstractRNG = Random.default_rng(), ) rps = reservoir_params(ps) - W = Float64.(rps.reservoir_matrix) - b = reservoir_bias(ps) - W_out = Float64.(readout_params(ps).weight) - C = resolve_conceptor(st, conceptor) - N = size(W, 1) - - x = init_state === nothing ? 0.5 .* randn(rng, N) : Float64.(init_state) - out_dims = size(W_out, 1) - Y = Matrix{Float64}(undef, out_dims, steps) - X = Matrix{Float64}(undef, N, steps) - - for n in 1:(washout + steps) - x = C * tanh.(W * x .+ b) - if n > washout - k = n - washout - @views X[:, k] = x - @views Y[:, k] = W_out * x + element_type = float(eltype(rps.reservoir_matrix)) + recurrent_weights = element_type.(rps.reservoir_matrix) + bias = reservoir_bias(ps) + readout_weights = element_type.(readout_params(ps).weight) + conceptor_matrix = element_type.(resolve_conceptor(st, conceptor)) + reservoir_dimension = size(recurrent_weights, 1) + checksquare(conceptor_matrix) + size(conceptor_matrix, 1) == reservoir_dimension || + throw( + DimensionMismatch( + "conceptor size must match the reservoir dimension $reservoir_dimension" + ) + ) + steps > 0 || throw(ArgumentError("steps must be positive")) + washout >= 0 || throw(ArgumentError("washout must be nonnegative")) + + state = init_state === nothing ? + element_type(0.5) .* randn(rng, element_type, reservoir_dimension) : + element_type.(init_state) + length(state) == reservoir_dimension || + throw(DimensionMismatch("init_state must have length $reservoir_dimension")) + output_dimension = size(readout_weights, 1) + outputs = zeros(element_type, output_dimension, steps) + states = zeros(element_type, reservoir_dimension, steps) + + for step in 1:(washout + steps) + state = conceptor_matrix * tanh.(recurrent_weights * state .+ bias) + if step > washout + output_index = step - washout + @views states[:, output_index] = state + @views outputs[:, output_index] = readout_weights * state end end - return Y, X + return outputs, states end -# ====================================================================== -# Conceptor morphing (Jaeger 2014, Section 3.2) -# ====================================================================== +# morphing @doc raw""" - morph_conceptor(st, weights) -> Matrix{Float64} + morph_conceptor(st, weights) -> Matrix Linearly combine stored conceptors into a morphed conceptor -``M = \sum_j \mu_j C_j`` (Jaeger 2014). `weights` is an iterable of `name => μ` -pairs, a `NamedTuple`, or a `Dict{Symbol,<:Real}`. Coefficients summing to one +(Jaeger 2014). `weights` is an iterable of name-to-weight pairs, a `NamedTuple`, +or a `Dict{Symbol,<:Real}`. Coefficients summing to one interpolate between the named prototypes; coefficients outside `[0, 1]` extrapolate. + +# Arguments + +- `st::NamedTuple`: State containing stored conceptors. +- `weights`: Named tuple, dictionary, or iterable of name-to-weight pairs. + +# Returns + +- The weighted sum of the named conceptors. + +# Throws + +- `KeyError`: If a named conceptor is absent. +- `ArgumentError`: If `weights` is empty. """ function morph_conceptor(st::NamedTuple, weights) - M = nothing - for (name, μ) in pairs(weights) - C = get_conceptor(st, Symbol(name)) - C === nothing && throw(KeyError(Symbol(name))) - contribution = Float64(μ) .* C - M = M === nothing ? contribution : M .+ contribution + morphed_conceptor = nothing + for (name, weight) in pairs(weights) + stored_conceptor = get_conceptor(st, Symbol(name)) + stored_conceptor === nothing && throw(KeyError(Symbol(name))) + contribution = convert(eltype(stored_conceptor), weight) .* stored_conceptor + morphed_conceptor = morphed_conceptor === nothing ? + contribution : morphed_conceptor .+ contribution end - M === nothing && throw(ArgumentError("no weights provided")) - return M + morphed_conceptor === nothing && throw(ArgumentError("no weights provided")) + return morphed_conceptor end # ====================================================================== @@ -569,8 +981,8 @@ end # ====================================================================== # Normalize a target into a row-major matrix (1 × T for a vector target). -_as_matrix(t::AbstractMatrix) = t -_as_matrix(t::AbstractVector) = reshape(t, 1, :) +_as_matrix(target::AbstractMatrix) = target +_as_matrix(target::AbstractVector) = reshape(target, 1, :) """ store_conceptors!(rng, concept, named_signals, ps, st; @@ -578,23 +990,39 @@ _as_matrix(t::AbstractVector) = reshape(t, 1, :) Derive and store one conceptor per named signal without modifying the reservoir weights. For each `name => signal` the reservoir is run from a cleared carry and -the conceptor `C(name) = R (R + α^{-2} I)^{-1}` of its state cloud is stored. Use +the conceptor of its state cloud is stored. Use this to build a conceptor library for classification, where each class pattern -gets its own conceptor. `aperture` is a scalar or a `name => α` dictionary. +gets its own conceptor. `aperture` is a scalar or a dictionary keyed by name. + +# Returns + +- The updated state. Reservoir parameters are unchanged. + +# See also + +[`load!`](@ref), [`train!`](@ref), [`store_conceptor!`](@ref) """ function store_conceptors!( - rng::AbstractRNG, concept::Conceptor, named_signals, ps::NamedTuple, st::NamedTuple; - aperture = 1.0, init_carry = nothing + rng::AbstractRNG, + concept::Conceptor, + named_signals, + ps::NamedTuple, + st::NamedTuple; + aperture = 1.0, + init_carry = nothing, ) st_acc = st for (name, signal) in named_signals - name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) sig = signal isa AbstractMatrix ? signal : reshape(signal, 1, :) st_reset = resetcarry!(rng, concept.model, st_acc.model; init_carry = init_carry) states, st_model = collectstates(concept.model, sig, ps.model, st_reset) st_acc = merge(st_acc, (; model = st_model)) - α = aperture_for(aperture, name) - store_conceptor!(st_acc, name, conceptor_from_states(states, α), α) + pattern_aperture = aperture_for(aperture, name) + store_conceptor!( + st_acc, name, conceptor_from_states(states, pattern_aperture), pattern_aperture + ) end return st_acc end @@ -610,18 +1038,48 @@ Train a single readout on conceptor-filtered reservoir features. For each are dropped, and features and the matching `named_targets` are pooled across patterns before the readout is fit with `train_method`. With `return_states=true` the pooled feature matrix is also returned. + +# Arguments + +- `named_signals`: Iterable of named input sequences. +- `named_targets`: Iterable containing one target sequence for every signal name. +- `train_method`: Readout fitting method; defaults to unregularized + [`StandardRidge`](@ref). + +# Keywords + +- `washout::Int = 0`: Initial columns removed from every feature and target. +- `return_states::Bool = false`: Return pooled conceptor-filtered features. +- `init_carry = nothing`: Initial carry passed when resetting each sequence. +- `kwargs...`: Additional options forwarded to `train`. + +# Returns + +- `(parameters, states)` by default. +- `((parameters, states), features)` when `return_states = true`. """ function train!( - rng::AbstractRNG, concept::Conceptor, named_signals, named_targets, - ps::NamedTuple, st::NamedTuple, train_method = StandardRidge(0.0); - washout::Int = 0, return_states::Bool = false, init_carry = nothing, kwargs... + rng::AbstractRNG, + concept::Conceptor, + named_signals, + named_targets, + ps::NamedTuple, + st::NamedTuple, + train_method = StandardRidge(0.0); + washout::Int = 0, + return_states::Bool = false, + init_carry = nothing, + kwargs..., + ) + targets = Dict{Symbol, AbstractMatrix}( + Symbol(name) => _as_matrix(target) for (name, target) in named_targets ) - targets = Dict{Symbol, AbstractMatrix}(Symbol(name) => _as_matrix(t) for (name, t) in named_targets) all_states = AbstractMatrix[] all_targets = AbstractMatrix[] for (name, signal) in named_signals - name isa Symbol || throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) haskey(targets, name) || throw(ArgumentError("no target data for pattern :$name")) st_reset = resetcarry!(rng, concept.model, st.model; init_carry = init_carry) st_tmp = merge(st, (; model = st_reset)) diff --git a/test/test_conceptors.jl b/test/test_conceptors.jl index f3fd426c..8b2abd92 100644 --- a/test/test_conceptors.jl +++ b/test/test_conceptors.jl @@ -5,68 +5,110 @@ using Statistics using ReservoirComputing # A small reproducible reservoir state cloud for algebra tests. -function sample_correlation(rng, N, L) - X = randn(rng, N, L) - return correlation_matrix(X) +function sample_correlation(rng, state_dimension, sample_count) + states = randn(rng, state_dimension, sample_count) + return correlation_matrix(states) end @testset "Conceptors" begin rng = Xoshiro(2024) @testset "conceptor matrix" begin - R = sample_correlation(rng, 20, 200) - C = conceptor_matrix(R, 10.0) - @test size(C) == (20, 20) - @test C ≈ C' # symmetric - ev = eigvals(Symmetric(C)) - @test all(>=(-1e-10), ev) # PSD - @test all(<=(1.0 + 1e-10), ev) # singular values ≤ 1 - @test conceptor_matrix(R, 1e8) ≈ I(20) atol = 1e-3 - @test norm(conceptor_matrix(R, 1e-6)) < 1e-6 - @test_throws ArgumentError conceptor_matrix(R, -1.0) + correlation = sample_correlation(rng, 20, 200) + conceptor = conceptor_matrix(correlation, 10.0) + @test size(conceptor) == (20, 20) + @test conceptor ≈ conceptor' + eigenvalues = eigvals(Symmetric(conceptor)) + @test all(>=(-1.0e-10), eigenvalues) + @test all(<=(1.0 + 1.0e-10), eigenvalues) + @test conceptor_matrix(correlation, 1.0e8) ≈ I(20) atol = 1.0e-3 + @test norm(conceptor_matrix(correlation, 1.0e-6)) < 1.0e-6 + @test_throws ArgumentError conceptor_matrix(correlation, -1.0) + @test_throws ArgumentError conceptor_matrix(correlation, Inf) + @test_throws ArgumentError correlation_matrix(zeros(3, 0)) + + for element_type in (Float32, Float64, BigFloat) + typed_states = element_type.(randn(rng, 6, 30)) + typed_correlation = correlation_matrix(typed_states) + typed_conceptor = conceptor_matrix(typed_correlation, 3) + @test eltype(typed_correlation) == element_type + @test eltype(typed_conceptor) == element_type + @test typed_conceptor ≈ typed_conceptor' + @test eltype(aperture_adapt(typed_conceptor, 2)) == element_type + end end @testset "aperture adaptation" begin - R = sample_correlation(rng, 15, 150) - C = conceptor_matrix(R, 5.0) - @test aperture_adapt(C, 1.0) ≈ C # γ = 1 is identity - # Prop 5: φ(φ(C, γ), β) = φ(C, γβ) - @test aperture_adapt(aperture_adapt(C, 2.0), 3.0) ≈ aperture_adapt(C, 6.0) - qs = [quota(aperture_adapt(C, g)) for g in (0.25, 1.0, 4.0, 16.0)] - @test issorted(qs) - P = aperture_adapt(C, Inf) - @test P * P ≈ P atol = 1e-8 # hardens to a projector - @test norm(aperture_adapt(C, 0.0)) < 1e-8 - @test conceptor_matrix(R, 7.0) ≈ aperture_adapt(conceptor_matrix(R, 1.0), 7.0) + correlation = sample_correlation(rng, 15, 150) + conceptor = conceptor_matrix(correlation, 5.0) + @test aperture_adapt(conceptor, 1.0) ≈ conceptor + # Repeated adaptations multiply their aperture factors. + @test aperture_adapt(aperture_adapt(conceptor, 2.0), 3.0) ≈ + aperture_adapt(conceptor, 6.0) + quotas = [ + quota(aperture_adapt(conceptor, factor)) for + factor in (0.25, 1.0, 4.0, 16.0) + ] + @test issorted(quotas) + projector = aperture_adapt(conceptor, Inf) + @test projector * projector ≈ projector atol = 1.0e-8 + @test norm(aperture_adapt(conceptor, 0.0)) < 1.0e-8 + @test conceptor_matrix(correlation, 7.0) ≈ + aperture_adapt(conceptor_matrix(correlation, 1.0), 7.0) end @testset "Boolean algebra" begin - C = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) - B = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) - A = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + first_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + second_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + third_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) - @test conceptor_not(conceptor_not(C)) ≈ C # double negation - @test conceptor_not(C) ≈ I(12) - C + @test conceptor_not(conceptor_not(first_conceptor)) ≈ first_conceptor + @test conceptor_not(first_conceptor) ≈ I(12) - first_conceptor # De Morgan - @test conceptor_or(C, B) ≈ conceptor_not(conceptor_and(conceptor_not(C), conceptor_not(B))) - @test conceptor_and(C, B) ≈ conceptor_not(conceptor_or(conceptor_not(C), conceptor_not(B))) + @test conceptor_or(first_conceptor, second_conceptor) ≈ + conceptor_not( + conceptor_and( + conceptor_not(first_conceptor), conceptor_not(second_conceptor) + ) + ) + @test conceptor_and(first_conceptor, second_conceptor) ≈ + conceptor_not( + conceptor_or( + conceptor_not(first_conceptor), conceptor_not(second_conceptor) + ) + ) # commutativity - @test conceptor_and(C, B) ≈ conceptor_and(B, C) - @test conceptor_or(C, B) ≈ conceptor_or(B, C) + @test conceptor_and(first_conceptor, second_conceptor) ≈ + conceptor_and(second_conceptor, first_conceptor) + @test conceptor_or(first_conceptor, second_conceptor) ≈ + conceptor_or(second_conceptor, first_conceptor) # AND shrinks, OR grows (quota ordering) - @test quota(conceptor_and(C, B)) <= quota(C) + 1e-8 - @test quota(C) <= quota(conceptor_or(C, B)) + 1e-8 + @test quota(conceptor_and(first_conceptor, second_conceptor)) <= + quota(first_conceptor) + 1.0e-8 + @test quota(first_conceptor) <= + quota(conceptor_or(first_conceptor, second_conceptor)) + 1.0e-8 # associativity (holds for conceptors, unlike idempotence) - @test conceptor_and(conceptor_and(C, B), A) ≈ conceptor_and(C, conceptor_and(B, A)) atol = 1e-6 - @test conceptor_or(conceptor_or(C, B), A) ≈ conceptor_or(C, conceptor_or(B, A)) atol = 1e-6 + @test conceptor_and( + conceptor_and(first_conceptor, second_conceptor), third_conceptor + ) ≈ conceptor_and( + first_conceptor, conceptor_and(second_conceptor, third_conceptor) + ) atol = 1.0e-6 + @test conceptor_or( + conceptor_or(first_conceptor, second_conceptor), third_conceptor + ) ≈ conceptor_or( + first_conceptor, conceptor_or(second_conceptor, third_conceptor) + ) atol = 1.0e-6 # neutral / absorbing elements - @test conceptor_and(C, Matrix(1.0I, 12, 12)) ≈ C atol = 1e-6 - @test conceptor_or(C, zeros(12, 12)) ≈ C atol = 1e-6 - # unexported infix sugar agrees - and_op, or_op, not_op = ReservoirComputing.:∧, ReservoirComputing.:∨, ReservoirComputing.:¬ - @test and_op(C, B) ≈ conceptor_and(C, B) - @test or_op(C, B) ≈ conceptor_or(C, B) - @test not_op(C) ≈ conceptor_not(C) + @test conceptor_and(first_conceptor, Matrix(1.0I, 12, 12)) ≈ + first_conceptor atol = 1.0e-6 + @test conceptor_or(first_conceptor, zeros(12, 12)) ≈ + first_conceptor atol = 1.0e-6 + + # The generalized definition must also work for singular conceptors. + first_singular_conceptor = Diagonal([0.8, 0.0, 0.4]) |> Matrix + second_singular_conceptor = Diagonal([0.5, 0.7, 0.0]) |> Matrix + @test conceptor_and(first_singular_conceptor, second_singular_conceptor) ≈ + Diagonal([1 / (1 / 0.8 + 1 / 0.5 - 1), 0, 0]) end @testset "wrapper / library" begin @@ -77,39 +119,55 @@ end st = initialstates(rng2, concept) @test haskey(ps, :model) @test isempty(st.conceptors) - C = conceptor_matrix(sample_correlation(rng2, 30, 100), 10.0) - store_conceptor!(st, :a, C, 10.0) - @test has_conceptor(st, :a) - @test get_conceptor(st, :a) ≈ C - st = set_active_conceptor(st, :a) - @test active_conceptor(st) ≈ C + conceptor = conceptor_matrix(sample_correlation(rng2, 30, 100), 10.0) + store_conceptor!(st, :first_pattern, conceptor, 10.0) + @test has_conceptor(st, :first_pattern) + @test get_conceptor(st, :first_pattern) ≈ conceptor + st = set_active_conceptor(st, :first_pattern) + @test active_conceptor(st) ≈ conceptor @test_throws KeyError set_active_conceptor(st, :missing) end @testset "load! / generate roundtrip" begin rng3 = Xoshiro(99) - bias_init(r, dims...) = 0.2f0 .* randn(r, Float32, dims...) - input_init(r, dims...) = 1.5f0 .* randn(r, Float32, dims...) - res_init(r, dims...) = rand_sparse(r, dims...; radius = 1.5f0, sparsity = 0.1f0) - esn = ESN(1, 80, 1; use_bias = true, init_bias = bias_init, - init_input = input_init, init_reservoir = res_init) + bias_init(random_generator, dims...) = + 0.2f0 .* randn(random_generator, Float32, dims...) + input_init(random_generator, dims...) = + 1.5f0 .* randn(random_generator, Float32, dims...) + res_init(random_generator, dims...) = + rand_sparse(random_generator, dims...; radius = 1.5f0, sparsity = 0.1f0) + esn = ESN( + 1, 80, 1; use_bias = true, init_bias = bias_init, + init_input = input_init, init_reservoir = res_init + ) concept = Conceptor(esn) ps = initialparameters(rng3, concept) st = initialstates(rng3, concept) - n = 1:1000 + sample_indices = 1:1000 period = 8.8342522 - sig = reshape(Float32.(sin.(2pi .* n ./ period)), 1, :) - ps, st = load!(rng3, concept, [:sine => sig], ps, st; aperture = 1000.0, washout = 400) + signal = reshape(Float32.(sin.(2pi .* sample_indices ./ period)), 1, :) + ps, st = load!( + rng3, concept, [:sine => signal], ps, st; aperture = 1000.0, washout = 400 + ) @test has_conceptor(st, :sine) - Y, X = generate(concept, ps, st; conceptor = :sine, steps = 200, washout = 300, rng = rng3) - @test size(Y) == (1, 200) - @test size(X) == (80, 200) - y = vec(Y) + outputs, states = generate( + concept, ps, st; conceptor = :sine, steps = 200, washout = 300, rng = rng3 + ) + @test size(outputs) == (1, 200) + @test size(states) == (80, 200) + output = vec(outputs) best = minimum( - sqrt(mean(abs2, y .- [s * sin(2pi * k / period + ph) for k in eachindex(y)])) / std(y) - for ph in 0:0.02:2pi, s in (-1.0, 1.0) + sqrt( + mean( + abs2, + output .- [ + sign * sin(2pi * index / period + phase) for + index in eachindex(output) + ], + ) + ) / std(output) for phase in 0:0.02:2pi, sign in (-1.0, 1.0) ) @test best < 0.05 @test count(>(0.5), conceptor_singular_values(get_conceptor(st, :sine))) < 80 @@ -120,12 +178,14 @@ end esn = ESN(1, 25, 1; use_bias = true) concept = Conceptor(esn) st = initialstates(rng4, concept) - Ca = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) - Cb = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) - store_conceptor!(st, :a, Ca, 10.0) - store_conceptor!(st, :b, Cb, 10.0) - @test morph_conceptor(st, (; a = 1.0, b = 0.0)) ≈ Ca - @test morph_conceptor(st, (; a = 0.3, b = 0.7)) ≈ 0.3 .* Ca .+ 0.7 .* Cb + first_conceptor = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + second_conceptor = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + store_conceptor!(st, :first_pattern, first_conceptor, 10.0) + store_conceptor!(st, :second_pattern, second_conceptor, 10.0) + @test morph_conceptor(st, (; first_pattern = 1.0, second_pattern = 0.0)) ≈ + first_conceptor + @test morph_conceptor(st, (; first_pattern = 0.3, second_pattern = 0.7)) ≈ + 0.3 .* first_conceptor .+ 0.7 .* second_conceptor @test_throws KeyError morph_conceptor(st, (; missing_name = 1.0)) end end