Commit e294c32f authored by Felix Kong's avatar Felix Kong

initial DMD class

parent 1dd79c0c
classdef dmd
%DMD Summary of this class goes here
% Detailed explanation goes here
properties (SetAccess = private)
X;
X2;
U0;
S0;
V0;
r; % rank to truncate to
U;
S;
V;
Atilde; % low-dimensional, truncated "weight-space" dynamics
W; % "weight-space" eigenvectors
Lambda;
Phi; % projection from weight-space to original data's space (e.g. image space)
verbose; % controls verbosity)
end
methods
function obj = dmd(data,varargin)
%DMD Construct an instance of this class
% The state at each timestep is a column vector)
obj.X = data(:,1:end-1);
obj.X2 = data(:,2:end);
obj.r = size(obj.X,2); % no truncation by default
if nargin == 1
obj.verbose = varargin{1};
else
obj.verbose = false;
end
obj.regression(); % it's fast, just run it on construction
end
function obj = regression(obj,varargin) % argument is truncation rank r
if nargin == 1
obj.r = varargin{1};
end
[obj.U0,obj.S0,obj.V0] = svd(obj.X,'econ');
% truncate
obj.U = obj.U0(:,1:obj_r);
obj.S = obj.S0(1:obj.r,1:obj.r);
obj.V = obj.V0(:,obj.r);
% construct weight-space dynamics
obj.Atilde = obj.U'*obj.X2*obj.V/(obj.S);
[obj.W,obj.Lambda] = eig(obj.Atilde);
% Phi projects from weight-space to original data space (e.g.
% image space)
obj.Phi = obj.X2*obj.V*inv(obj.S)*obj.W;
end
function X2hat = reconstruct(obj)
% reconstructs an estimate of X2 using the low-rank
% "weight-space" dynamics x_r(k+1) = Atilde*x_r(k)
%
% useful for getting a feel for how accurate the low-rank
% approximation is.
% UNTESTED, NOT SURE IF THIS IS MATHEMATICALLY CORRECT
X2hat = obj.Phi*obj.Lambda*obj.Phi\obj.X;
end
function xhat = predict(obj,initial_condition,num_timesteps)
assert(size(initial_condition) == size(obj.X,1));
xhat = nan(length(initial_condition),num_timesteps);
xhat(:,1) = obj.Phi\initial_condition;
for k = 2:num_timesteps
xhat(:,k) = obj.Phi*obj.Lambda*obj.Phi\xhat(:,k-1);
end
end
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment