Paper discussion

Date: November 03, 2024

These are rough lecture notes for a paper review discussion held spring 2024 with the SMU data Applied Mathematics/Data Science Journal Club. This work is simply a partial account of the content in the Cited paper.

Author: Mason A. McCallum

Citation

Bakarji J, Champion K, Nathan Kutz J, Brunton SL. 2023 Discovering governing equations from partial measurements with deep delay autoencoders. Proc. R. Soc. A 479: 20230422. https://doi.org/10.1098/rspa.2023.0422

Introduction to SINDY

Here we will empirically determine a mathematical model of observed dynamical variables. To employ this technique we must collect time series data and coorisponding time derivative data. Let $x:\mathbb{R}^+ \rightarrow \mathbb{R}^d$

$X=[x(t_1), x(t_2),...,x(t_m)]^T$ $\dot{X}=[\dot{x}(t_1), \dot{x}(t_2),...,\dot{x}(t_m)]^T$

These times series matrices $X$ and $\dot{X}$ are $\mathbb{R}^{m\times d}$

We will regularize our optimization to promote models with as few terms as possible.

The continous time dynamical system is written as $\frac{d}{dt}x(t)=f(x)$ the discritized equivelant of this system via our time series data is $\dot{X}=\Theta(X)\xi$

Definition: The Library $\Theta$ is a large matrix $\Theta(X)\in\mathbb{R}^{m\times s}$ $\Theta(X)=[1\;\; X_1\;\;X_2\;\;X_3\;\; X_1X_2\;\; X_1X_3\;\; X_3X_2\;\; X_1^2\cdot\cdot\cdot X^d \cdot\cdot\cdot sin(X)\cdot\cdot\cdot ]$ Each column in the above definition corresponds to applying some mapping to our time series data.

Definition: The coefincent matrix $\xi \in \mathbb{R}^{s \times d}$ is a sparce matrix which we will optimize to satisfy $\dot{X}=\Theta(X)\xi$. This matrix serves as a means of selecting which library functions belong in the model.

\[L_{SINDY}(\xi')=||\dot{X}-\Theta(X)\xi'||_2+\lambda||\xi'||_1\]

using Plots
using LinearAlgebra

FT = Float64

# define the Rossler attractor
@kwdef mutable struct Rossler
    dt::FT = 0.02
    a::FT = 0.2
    b::FT = 0.2
    c::FT = 5.7
    x::FT = 0.1
    y::FT = 1.0
    z::FT = 0.0
    dx::FT = 0.0
    dy::FT = 0.0
    dz::FT = 0.0
end

function step!(R::Rossler)
    R.dx = -R.y - R.z;             
    R.dy = R.x + R.a*R.y;       
    R.dz = R.b + R.z*(R.x-R.c);     

    R.x += R.dt * R.dx
    R.y += R.dt * R.dy
    R.z += R.dt * R.dz
end

function keytostring(keyline)
  numVars = size(keyline)[1]
  for var in 1:numVars
    print("*x",var,"^",keyline[var])
  end
end 

function printSystem(key,coefs)
  sizeLib, numVars = size(coefs)
  for var in 1:numVars
    print("\n ẋ",var,"=")
    for libterm in 1:sizeLib
      if !(coefs[libterm, var] ≈ 0.0)
        print("+ ",coefs[libterm, var],"(")
        keytostring(key[libterm,:])
        print(")")
      end
    end
  end
end

printSystem (generic function with 1 method)

Build data matrices

numSteps = 10000
attractor = Rossler()

# Time series data collection
X = zeros(FT, numSteps, 3) #ℜᵐˣ³
Ẋ = zeros(FT, numSteps, 3) #ℜᵐˣ³

for ti=1:numSteps
    X[ti, :] = Vector{FT}([attractor.x, attractor.y, attractor.z])
    Ẋ[ti, :] = Vector{FT}([attractor.dx, attractor.dy, attractor.dz])
    step!(attractor)
end

# Library Formation
funcs = [(x,y,z)->x^a*y^b*z^c for c=0:2 for b=0:2 for a=0:2];
labels = [(a,b,c) for c=0:2 for b=0:2 for a=0:2]
order = sortperm(labels, by = x->sum(x)*sum(x.!=0))
funcs = funcs[order]
labels = labels[order]

Θ = zeros(FT, numSteps, length(funcs));
xi = zeros(FT, length(funcs), 3)# Ξ

for (idx, func) in enumerate(funcs)
  Θ[:,idx] = func.(X[:,1], X[:,2], X[:,3])
end

ptrue = plot(X[:,1],X[:,2],X[:,3],legend=false, title="truth");

Optimize

$\begin{align*} &\dot{x}=-y-z\\ &\dot{y}=x+0.2y\\ &\dot{z}=0.2+z(x-5.7)\\ \end{align*}$

"""
SparsifyDynamics for finding ξ
Theta: Θ(X)
dXdt: Ẋ
lambda: Regularization parameter weighting sparsinty of ξ. 
"""
function sparsifyDynamics(Theta::Matrix{FT}, dXdt::Matrix{FT}, lambda::FT)
  Xi = Theta\dXdt   # Initial guess: Least-squares
  
  for k=1:10
    smallinds = (abs.(Xi) .< lambda)
    Xi[smallinds] .= 0  # Threasholding
    for col = 1:size(dXdt,2)
      biginds = .!smallinds[:,col]
      Xi[biginds,col] = Theta[:,biginds]\dXdt[:,col]
    end
  end
  return Xi
end

# Now we seek the optimization of 
# Ξ = argmin||Ẋ-Θ(X)Ξ||₂ + λ||Ξ||₁
λ = 0.055
xi = sparsifyDynamics(Θ, Ẋ, λ)  # Here one could use Lasso regression among other sparse optimizers

# Now to print out the system
coefs = xi
key = transpose(reinterpret(reshape, Int, labels))
printSystem(key,coefs)

 ẋ1=+ -1.0054302319769184(*x1^0*x2^1*x3^0)+ -0.9871832765872768(*x1^0*x2^0*x3^1)
 ẋ2=+ 0.9998938754508867(*x1^1*x2^0*x3^0)+ 0.22139086315654422(*x1^0*x2^1*x3^0)
 ẋ3=+ -5.765816978762432(*x1^0*x2^0*x3^1)+ 1.0220655880065033(*x1^1*x2^0*x3^1)+ 0.08193878455408737(*x1^0*x2^1*x3^1)

function fmodel(x,y,z)
    dX = [0.0, 0.0, 0.0]
    for col in 1:size(xi,2)
        for row in 1:size(xi,1)
            dX[col] += xi[row,col]*funcs[row](x,y,z)
        end
    end
    return dX
end

xmodel = zeros(FT, numSteps, 3)
xmodel[1,:] = [0.01, 1.0, 0.0]
for i=2:numSteps
    xmodel[i,:] = xmodel[i-1,:] +  0.02 * fmodel(xmodel[i-1,1], xmodel[i-1,2], xmodel[i-1,3])
end
pmodel = plot(xmodel[:,1],xmodel[:,2],xmodel[:,3],legend=false, title="model");
plot(ptrue, pmodel, layout=(1,2))

<?xml version=”1.0” encoding=”utf-8”?>

In the above discussion we collected data in the coordinate system of the dynamical variables x,y,z. Suppose however the observed data were not in the dynamical coordinate system. This is akin to developing orbital mechanics before Copernicus suggested a heliocentric coordinate system. The differential equations based on a geocentric coordinate system are very messy. Once Copernicus posited a heliocentric coordinate system the dynamical system became discoverable and led to Kepler’s findings. In the following discussion we will observe data from a system. Then in the spirit of Copernicus discover an appropriate coordinate system for the dynamics this will be accomplished using an autoencoder to map our observables into a latent space. Then within that latent space SINDy will play the role of Kepler to discover a dynamical model. Aside from the very general conceptual detail of mathematical modeling in an appropriate coordinate system the interesting details of the autoencoder is the selection of loss functions to constrain the optimization and discover the appropriate mapping from observation space into the dynamical space.

consider the nonlinear differential equation $\frac{d}{dt}x(t)=f(x(t))$ Let $y\in\mathbb{R}^d$ be noisy measurments of the system via some observation operator. $y(t)=g(x(t))+\eta$. For example let g be the projection of the system onto $e_1$. Consider some unkown mapping from measure time delay space to the dynamical variables $z(t)=\phi(y(t))$. Appon finding $\phi$ (via an autoencoder) the goal is to discover the dynamical system $\frac{d}{dt}z(t)=h(z(t))$ (SINDy)

Time series data will be collected into the Hankel matrix $H= \begin{bmatrix} &y(t_1) &y(t_2) &\cdot\cdot\cdot &y(t_q)\\ &y(t_2) &y(t_3) &\cdot\cdot\cdot &y(t_{q+1})\\ &\cdot\cdot\cdot &\cdot\cdot\cdot &\cdot\cdot\cdot &\cdot\cdot\cdot\\ &y(t_n) &y(t_{n+1}) &\cdot\cdot\cdot &y(t_{n+q+1})\\ \end{bmatrix}=[h_1 h_2 \cdot\cdot\cdot h_q]$

The autoencoder network is a mapping from the column space of the H matrix into our latent space. The punchline is from noisy observation data (e.g. just collecting a single dynamical variable) the full dynamics can be discovered after finding the appropriate mapping into latent space. Note that in the paper one can account for noisy data by taking the reduced SVD of the H matrix.

Warning

This code is not finished. Not all the Loss functions are implemented therefore our search for the appropriate mapping to latent space is not going to be correct. Most of the bones are there but the loss functions are needed.

using Flux, Plots, ProgressMeter

# Time series data collection
X_sim = zeros(FT, numSteps, 3)
Ẋ_sim = zeros(FT, numSteps, 3)
for ti=1:numSteps
    step!(attractor)
    X_sim[ti, :] = Vector{FT}([attractor.x, attractor.y, attractor.z])
    Ẋ_sim[ti, :] = Vector{FT}([attractor.dx, attractor.dy, attractor.dz])
end

function hankel(data, d)
  len = length(data)
  q = len - d
  H = Matrix{typeof(data[1])}(undef, d, q)
  for col=1:q
    H[:, col] = data[col:col+d-1]
  end
  return H
end

n = 60
n2 = 30
m = 15
r = 3
H = hankel(X_sim[:,1], n)
data = [Vector(H[:,col]) for col=1:size(H)[2]]
enc = Chain(
            Dense(n => n2),
            Dense(n2 => m), 
            Dense(m => m), 
            Dense(m => m), 
            Dense(m => m), 
            Dense(m => r), 
          )
dec = Chain(
            Dense(r => m),
            Dense(m => m), 
            Dense(m => m), 
            Dense(m => m), 
            Dense(m => n2), 
            Dense(n2 => n)
          )
model = Chain(enc, dec)

optim = Flux.setup(Adam(), model)

@showprogress for epoch=1:100
  Flux.train!(model, data[1:5000], optim) do m, x
    sum((m(x).-x).^2)+(m[1](x)[1]-x[1])^2
  end
end

[32mProgress: 100%|█████████████████████████████████████████| Time: 0:01:12[39m

latent_data = [model[1](d) for d in data]
X = latent_data
Ẋ = [(latent_data[idx]-latent_data[idx+1])/0.02 for idx=1:length(X)-1]
X = transpose(hcat(X[1:end-1]...))
Ẋ = Matrix(transpose(hcat(Ẋ...)))

attractor = Rossler()

# Library Formation
funcs = [(x,y,z)->x^a*y^b*z^c for c=0:2 for b=0:2 for a=0:2];
labels = [(a,b,c) for c=0:2 for b=0:2 for a=0:2]
order = sortperm(labels, by = x->sum(x)*sum(x.!=0))
funcs = funcs[order]
labels = labels[order]
numSteps = size(X)[1]
Θ = zeros(FT, numSteps, length(funcs));
xi = zeros(FT, length(funcs), 3)# Ξ

for (idx, func) in enumerate(funcs)
  Θ[:,idx] = func.(X[:,1], X[:,2], X[:,3])
end

# Now we seek the optimization of 
# Ξ = argmin||Ẋ-Θ(X)Ξ||₂ + λ||Ξ||₁
λ = 0.061
xi = sparsifyDynamics(Θ, Ẋ, λ)

coefs = xi
key = transpose(reinterpret(reshape, Int, labels));
printSystem(key,coefs)

function fmodel(x,y,z)
    dX = [0.0, 0.0, 0.0]
    for col in 1:size(xi,2)
        for row in 1:size(xi,1)
            dX[col] += xi[row,col]*funcs[row](x,y,z)
        end
    end
    return dX
end

xmodel = zeros(FT, numSteps, 3)
xmodel[1,:] = [0.01, 1.0, 0.0]
for i=2:numSteps
    xmodel[i,:] = xmodel[i-1,:] +  0.02 * fmodel(xmodel[i-1,1], xmodel[i-1,2], xmodel[i-1,3])
end
pmodel = plot(xmodel[:,1],xmodel[:,2],xmodel[:,3],legend=false, title="model");
plot(ptrue, pmodel, layout=(1,2))

 ẋ1=+ 0.6963388095479174(*x1^0*x2^0*x3^0)+ 0.2792826819737061(*x1^1*x2^0*x3^0)+ 0.061130985589659785(*x1^0*x2^0*x3^1)
 ẋ2=+ 39.65671143525014(*x1^0*x2^0*x3^0)+ -2.056409520378554(*x1^1*x2^0*x3^0)+ 1.8252778746922187(*x1^0*x2^1*x3^0)+ 5.282539361932167(*x1^0*x2^0*x3^1)+ 4.255475137871766(*x1^2*x2^0*x3^0)+ -0.1829295578441074(*x1^1*x2^1*x3^0)+ -0.49365590871893644(*x1^1*x2^0*x3^1)
 ẋ3=+ -115.7144459596724(*x1^0*x2^0*x3^0)+ 135.70854614965396(*x1^1*x2^0*x3^0)+ -5.519357602826046(*x1^0*x2^1*x3^0)+ -10.005406528755964(*x1^0*x2^0*x3^1)+ -16.16325533002859(*x1^2*x2^0*x3^0)+ 0.6774592686814631(*x1^1*x2^1*x3^0)+ 1.307410828751469(*x1^1*x2^0*x3^1)