include online DMD, window DMD

README.md

@@ -1,2 +1,15 @@
 # dmdtools
 A library of tools for computing variants of Dynamic Mode Decomposition
+
+## algorithms implemented in Matlab
+1. Standard batch processed total-least-squares DMD (TDMD)  
+2. Streaming DMD (including total-least-squares)  
+3. Online DMD (including weighted version)  
+4. Window DMD (including weighted version)    
+
+## algorithms implemented in Python
+1. Standard batch processed DMD  
+2. Kernel DMD with polynomial kernel  
+3. Streaming DMD  
+4. Online DMD (including weighted version)  
+5. Window DMD (including weighted version)  

matlab/OnlineDMD.m

@@ -0,0 +1,130 @@
+% OnlineDMD is a class that implements online dynamic mode decomposition
+% The time complexity (multiply-add operation for one iteration) is O(4n^2)
+% , and space complexity is O(2n^2), where n is the state dimension.
+%
+% Algorithm description:
+%       At time step k, define two matrix X(k) = [x(1),x(2),...,x(k)], 
+%       Y(k) = [y(1),y(2),...,y(k)], that contain all the past snapshot 
+%       pairs, where x(k), y(k) are the n dimensional state vector, 
+%       y(k) = f(x(k)) is the image of x(k), f() is the dynamics. 
+%       Here, if the (discrete-time) dynamics are given by z(k) = f(z(k-1))
+%       , then x(k), y(k) should be measurements corresponding to 
+%       consecutive states z(k-1) and z(k).
+%       At time step k+1, we need to include new snapshot pair x(k+1), y(k+1)
+%       We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively 
+%       by efficient rank-1 updating online DMD algorithm.
+%       An exponential weighting factor can be used to place more weight on
+%       recent data.
+%
+% Usage:
+%       odmd = OnlineDMD(n,weighting)
+%       odmd.initialize(Xq,Yq)
+%       odmd.initilizeghost()
+%       odmd.update(x,y)
+%       [evals, modes] = odmd.computemodes()
+%        
+% properties:
+%       n: state dimension
+%       weighting: weighting factor in (0,1]
+%       timestep: number of snapshot pairs processed
+%       A: DMD matrix, size n by n
+%       P: matrix that contains information about past snapshots, size n by n
+% 
+% methods:
+%       initialize(Xq, Yq), initialize online DMD algorithm with q snapshot
+%                           pairs stored in (Xq, Yq)
+%       initializeghost(), initialize online DMD algorithm with epsilon 
+%                           small (1e-15) ghost snapshot pairs before t=0
+%       update(x,y), update when new snapshot pair (x,y) becomes available
+%                   Here, if the (discrete-time) dynamics are given by 
+%                   z(k) = f(z(k-1)), then (x,y) should be measurements 
+%                   correponding to consecutive states z(k-1) and z(k).
+%       computemodes(), compute and return DMD eigenvalues and DMD modes
+%
+% Authors: 
+%   Hao Zhang
+%   Clarence W. Rowley
+% 
+% Reference:
+% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta,
+% "Online Dynamic Mode Decomposition for Time-varying Systems,"
+% arXiv preprint arXiv:1707.02876, 2017.
+% 
+% Created:
+%   April 2017.
+%
+% To look up the documentation in the command window, type help OnlineDMD
+
+classdef OnlineDMD < handle
+    properties
+        n = 0;                      % state dimension
+        weighting = 1;                 % weighting factor in (0,1]
+        timestep = 0;               % number of snapshots processed
+        A;          % DMD matrix
+        P;          % matrix that contains information about past snapshots
+    end
+
+    methods
+        function obj = OnlineDMD(n,weighting)
+            % Creat an object for online DMD
+            % Usage: odmd = OnlineDMD(n,weighting)
+            if nargin == 2
+                obj.n = n;
+                obj.weighting = weighting;
+                obj.timestep = 0;
+                obj.A = zeros(n,n);
+                obj.P = zeros(n,n);
+            end
+        end
+
+        function initialize(obj, Xq, Yq)
+            % Initialize OnlineDMD with q snapshot pairs stored in (Xq, Yq)
+            % Usage: odmd.initialize(Xq,Yq)
+            q = length(Xq(1,:));
+            if(obj.timestep == 0 && rank(Xq) == obj.n)
+                weight = (sqrt(obj.weighting)).^(q-1:-1:0);
+                Xq = Xq.*weight;
+                Yq = Yq.*weight;
+                obj.A = Yq*pinv(Xq);
+                obj.P = inv(Xq*Xq')/obj.weighting;
+            end
+            obj.timestep = obj.timestep + q;
+        end
+
+        function initializeghost(obj)
+            % Initialize online DMD with epsilon small (1e-15) ghost 
+            % snapshot pairs before t=0
+            % Usage: odmd.initilizeghost()
+            epsilon = 1e-15;
+            obj.A = randn(obj.n, obj.n);
+            obj.P = (1/epsilon)*eye(obj.n);
+        end
+
+        function update(obj, x, y)
+            % Update the DMD computation with a new pair of snapshots (x,y)
+            % Here, if the (discrete-time) dynamics are given by z(k) = 
+            % f(z(k-1)), then (x,y) should be measurements correponding to 
+            % consecutive states z(k-1) and z(k).
+            % Usage: odmd.update(x, y)
+
+            % compute P*x matrix vector product beforehand
+            Px = obj.P*x;
+            % Compute gamma
+            gamma = 1/(1+x'*Px);
+            % Update A
+            obj.A = obj.A + (gamma*(y-obj.A*x))*Px';
+            % Update P, group Px*Px' to ensure positive definite
+            obj.P = (obj.P - gamma*(Px*Px'))/obj.weighting;
+            % ensure P is SPD by taking its symmetric part
+            obj.P = (obj.P+(obj.P)')/2;
+            % time step + 1
+            obj.timestep = obj.timestep + 1;
+        end
+
+        function [evals, modes] = computemodes(obj)
+            % Compute DMD eigenvalues and DMD modes at current time
+            % Usage: [evals, modes] = odmd.modes()
+            [modes, evals] = eig(obj.A, 'vector');
+        end
+    end
+end

matlab/README.md

@@ -0,0 +1,15 @@
+# README for Matlab implementation of variants of Dynamic Mode Decomposition
+
+## Implementation
+1.**tdmd.m** implements the total-least-squares DMD function.  
+2.**StreamingDMD.m** implements **StramingDMD** class.  
+3.**StreamingTDMD.m** implements **StreamingTDMD** class.  
+4.**OnlineDMD.m** implements **OnlineDMD** class.  
+5.**WindomDMD.m** implements **WindowDMD** class.  
+
+## Demos
+1.**tdmd_run.m** demostrates the use of total-least-squares DMD function.  
+2.**sdmd_run.m** demostrates the use of **StreamingDMD** class.  
+3.**stdmd_run.m** demostrates the use of **StreamingTDMD** class.  
+4.**online_demo.m** demostrates the use of **OnlineDMD** class.  
+5.**window_demo.m** demostrates the use of **WindowDMD** class.  

matlab/WindowDMD.m

@@ -0,0 +1,156 @@
+% WindowDMD is a class that implements window dynamic mode decomposition
+% The time complexity (multiply-add operation for one iteration) is
+% O(8n^2), and space complexity is O(2wn+2n^2), where n is the state
+% dimension, w is the window size.
+% 
+% Algorithm description:
+%       At time step k, define two matrix 
+%       X(k) = [x(k-w+1),x(k-w+2),...,x(k)], Y(k) = [y(k-w+1),y(k-w+2),
+%       ...,y(k)], that contain the recent w snapshot pairs from a finite 
+%       time window, where x(k), y(k) are the n dimensional state vector, 
+%       y(k) = f(x(k)) is the image of x(k), f() is the dynamics. 
+%       Here, if the (discrete-time) dynamics are given by z(k) = f(z(k-1))
+%       , then x(k), y(k) should be measurements corresponding to 
+%       consecutive states z(k-1) and z(k).
+%       At time step k+1, we need to forget old snapshot pair xold = 
+%       x(k-w+1), yold = y(k-w+1), and remember new snapshot pair xnew = 
+%       x(k+1), ynew = y(k+1).
+%       We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively 
+%       by efficient rank-2 updating window DMD algorithm.
+%       An exponential weighting factor can be used to place more weight on
+%       recent data.
+%        
+% Usage:
+%       wdmd = WindowDMD(n,w,weighting)
+%       wdmd.initialize(Xw,Yw)
+%       wdmd.update(xnew, ynew)
+%       [evals, modes] = wdmd.computemodes()
+%
+% properties:
+%       n: state dimension
+%       w: finite time window size
+%       weighting: weighting factor in (0,1]
+%       timestep: number of snapshot pairs processed
+%       Xw: recent w snapshots x stored in Xw, size n by w
+%       Yw: recent w snapshots y stored in Yw, size n by w
+%       A: DMD matrix for w snapshot pairs, size n by n
+%       P: Matrix that contains information about recent w snapshots, 
+%          size n by n
+% 
+% methods:
+%       initialize(Xw, Yw), initialize window DMD algorithm
+%       update(xnew, ynew), update by forgetting old snapshot pairs, 
+%       and remeber new snapshot pair, move sliding window forward
+%       At time k+1, X(k+1) = [x(k-w+2),x(k-w+2),...,x(k+1)], 
+%       Y(k+1) = [y(k-w+2),y(k-w+2),...,y(k+1)], 
+%       we should take xnew = x(k+1), ynew = y(k+1)
+%       computemodes(), compute and return DMD eigenvalues and DMD modes
+%
+% Authors: 
+%   Hao Zhang
+%   Clarence W. Rowley
+% 
+% Reference:
+% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta,
+% "Online Dynamic Mode Decomposition for Time-varying Systems,"
+% arXiv preprint arXiv:1707.02876, 2017.
+% 
+% Created:
+%   April 2017.
+%
+% To look up the documentation, type help WindowDMD
+
+classdef WindowDMD < handle
+    properties
+        n = 0;              % state dimension
+        w = 0;              % window size
+        weighting = 1;      % weighting factor in (0,1]
+        timestep = 0;       % number of snapshots processed
+        Xw;     % recent w snapshots x stored in matrix Xw
+        Yw;     % recent w snapshots y stored in matrix Yw
+        A;      % DMD matrix for w snapshot pairs, size n by n
+        P;      % Matrix that contains information about recent w snapshots
+                % , size n by n
+    end
+
+    methods
+        function obj = WindowDMD(n, w, weighting)
+            % Creat an object for window DMD
+            % Usage: wdmd = WindowDMD(n,w,weighting)
+            if nargin == 3
+                obj.n = n;
+                obj.w = w;
+                obj.weighting = weighting;
+                obj.timestep = 0;
+                obj.Xw = zeros(n,w);
+                obj.Yw = zeros(n,w);
+                obj.A = zeros(n,n);
+                obj.P = zeros(n,n);
+            end
+        end
+
+        function initialize(obj, Xw, Yw)
+            % Initialize WindowDMD with w snapshot pairs stored in (Xw, Yw)
+            % Usage: wdmd.initialize(Xw,Yw)
+
+            % initialize Xw, Yw
+            obj.Xw = Xw; obj.Yw = Yw;
+            % initialize A, P
+            q = length(Xw(1,:));
+            if(obj.timestep == 0 && obj.w == q && rank(Xw) == obj.n)
+                weight = (sqrt(obj.weighting)).^(q-1:-1:0);
+                Xw = Xw.*weight;
+                Yw = Yw.*weight;
+                obj.A = Yw*pinv(Xw);
+                obj.P = inv(Xw*Xw')/obj.weighting;
+            end
+            obj.timestep = obj.timestep + q;
+        end
+
+        function update(obj, xnew, ynew)    
+            % Update the DMD computation by sliding the finite time window 
+            % forward.
+            % Forget the oldest pair of snapshots (xold, yold), and 
+            % remembers the newest pair of snapshots (xnew, ynew) in the 
+            % new time window. If the new finite time window at time step 
+            % k+1 includes recent w snapshot pairs as
+            % X(k+1) = [x(k-w+2),x(k-w+3),...,x(k+1)], 
+            % Y(k+1) = [y(k-w+2),y(k-w+3),...,y(k+1)],
+            % where y(k) = f(x(k)) and f is the dynamics, then we should 
+            % take xnew = x(k+1), ynew = y(k+1)
+            % Usage: wdmd.update(xnew, ynew)
+
+            % define old snapshots to be discarded
+            xold = obj.Xw(:,1); yold = obj.Yw(:,1);
+            % Update recent w snapshots
+            obj.Xw = [obj.Xw(:,2:end), xnew];
+            obj.Yw = [obj.Yw(:,2:end), ynew];
+
+            % direct rank-2 update
+            % define matrices
+            U = [xold, xnew]; V = [yold, ynew]; 
+            C = diag([-(obj.weighting)^(obj.w),1]);
+            % compute PkU matrix matrix product beforehand
+            PkU = obj.P*U;
+            % compute AkU matrix matrix product beforehand
+            AkU = obj.A*U;
+            % compute Gamma
+            Gamma = inv(inv(C)+U'*PkU);
+            % update A
+            obj.A = obj.A + (V-AkU)*(Gamma*PkU');
+            % update P
+            obj.P = (obj.P - PkU*(Gamma*PkU'))/obj.weighting;
+            % ensure P is SPD by taking its symmetric part
+            obj.P = (obj.P+(obj.P)')/2;
+
+            % time step + 1
+            obj.timestep = obj.timestep + 1;
+        end
+
+        function [evals, modes] = computemodes(obj)
+            % Compute DMD eigenvalues and DMD modes at current time step
+            % Usage: [evals, modes] = wdmd.computemodes()
+            [modes, evals] = eig(obj.A, 'vector');
+        end
+    end
+end