From 6173b32c3780a2986a417ccb04950f86463db5eb Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 16:59:12 -0400 Subject: [PATCH 01/32] add online DMD and window DMD --- matlab/OnlineDMD.m | 114 ++++++++++++++++++++++++++++++ matlab/WindowDMD.m | 113 ++++++++++++++++++++++++++++++ matlab/online_demo.m | 117 +++++++++++++++++++++++++++++++ matlab/window_demo.m | 103 +++++++++++++++++++++++++++ python/dmdtools/online.py | 109 +++++++++++++++++++++++++++++ python/dmdtools/window.py | 112 +++++++++++++++++++++++++++++ python/scripts/online_demo.py | 128 ++++++++++++++++++++++++++++++++++ python/scripts/window_demo.py | 114 ++++++++++++++++++++++++++++++ 8 files changed, 910 insertions(+) create mode 100644 matlab/OnlineDMD.m create mode 100644 matlab/WindowDMD.m create mode 100644 matlab/online_demo.m create mode 100644 matlab/window_demo.m create mode 100644 python/dmdtools/online.py create mode 100644 python/dmdtools/window.py create mode 100644 python/scripts/online_demo.py create mode 100644 python/scripts/window_demo.py diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m new file mode 100644 index 0000000..a15772b --- /dev/null +++ b/matlab/OnlineDMD.m @@ -0,0 +1,114 @@ +% OnlineDMD is a class that implements online dynamic mode decomposition +% The time complexity (for one iteration) is O(n^2), and space complexity is +% O(n^2), where n is the state dimension. +% +% Algorithm description: +% At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-1 updating online DMD algorithm. +% +% Usage: +% odmd = OnlineDMD(n,lambda) +% odmd.initialize(Xq,Yq) +% odmd.initilizeghost() +% odmd.update(x,y) +% [evals, modes] = odmd.computemodes() +% +% properties: +% n: state dimension +% lambda: weighting factor between 0 and 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 mdoes +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Created: +% April 2017. +% +% To look up the documentation in the command window, type help OnlineDMD + +classdef OnlineDMD < handle + properties + n = 0; % state dimension + lambda = 1; % weighting factor + timestep = 0; % number of snapshots processed + A; % DMD matrix + P; % matrix that contains information about past snapshots + end + + methods + function obj = OnlineDMD(n,lambda) + % Creat an object for online DMD + % Usage: odmd = OnlineDMD(n,lambda) + if nargin == 2 + obj.n = n; + obj.lambda = lambda; + 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 && q>=obj.n) + sqrtlambda = sqrt(obj.lambda); + for i = 1:q + Xq(:,i) = Xq(:,i)*sqrtlambda^(q-i); + Yq(:,i) = Yq(:,i)*sqrtlambda^(q-i); + end + obj.A = Yq*pinv(Xq); + obj.P = inv(Xq*Xq')/obj.lambda; + 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; + alpha = 1.0/epsilon; + obj.A = randn(obj.n, obj.n); + obj.P = alpha*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 gamma + gamma = 1/(1+x'*(obj.P*x)); + % Update A + obj.A = obj.A + gamma*((y-obj.A*x)*(x'*obj.P)); + % Update P + obj.P = (obj.P - gamma*((obj.P*x)*(x'*obj.P)))/obj.lambda; + + obj.timestep = obj.timestep + 1; + end + + function [evals, modes] = computemodes(obj) + % Compute and return DMD eigenvalues and DMD modes at current time + % Usage: [evals, modes] = odmd.modes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m new file mode 100644 index 0000000..b13c143 --- /dev/null +++ b/matlab/WindowDMD.m @@ -0,0 +1,113 @@ +% WindowDMD is a class that implements window dynamic mode decomposition +% The time complexity (for one iteration) is O(n^2), and space complexity is +% O(n^2), where n is the state dimension +% +% Algorithm description: +% At time step k, define two matrix +% Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-2 updating window DMD algorithm. +% +% Usage: +% wdmd = WindowDMD(n,w) +% wdmd.initialize(Xq,Yq) +% wdmd.update(xold, yold, xnew, ynew) +% [evals, modes] = wdmd.computemodes() +% +% properties: +% n: state dimension +% w: finite time window size +% timestep: number of snapshot pairs processed +% A: Intermediate DMD matrix for w-1 snapshot pairs, size n by n +% B: Intermediate DMD matrix for w-1 snapshot pairs, size n by n +% M: Matrix that contains information about recent w-1 snapshots, size n by n +% +% methods: +% initialize(Xq, Yq), initialize window DMD algorithm +% update(xold, yold, xnew, ynew), update when new snapshot pair becomes available +% computemodes(), compute and return DMD eigenvalues and DMD mdoes +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Created: +% April 2017. +% +% To look up the documentation, type help WindowDMD + +classdef WindowDMD < handle + properties + n = 0; % state dimension + w = 0; % weighting factor + timestep = 0; % number of snapshots processed + A; % Intermediate DMD matrix for w-1 snapshot pairs, size n by n + B; % Intermediate DMD matrix for w-1 snapshot pairs, size n by n + M; % Matrix that contains information about recent w-1 snapshots, size n by n + end + + methods + function obj = WindowDMD(n,w) + % Creat an object for window DMD + % Usage: wdmd = WindowDMD(n,w) + if nargin == 2 + obj.n = n; + obj.w = w; + obj.A = zeros(n,n); + obj.B = zeros(n,n); + obj.M = zeros(n,n); + end + end + + function initialize(obj, Xq, Yq) + % Initialize WnlineDMD with q snapshot pairs stored in (Xq, Yq) + % Usage: wdmd.initialize(Xq,Yq) + q = length(Xq(1,:)); + if(obj.timestep == 0 && obj.w == q && obj.w >= obj.n+1) + obj.A = Yq*pinv(Xq); + obj.B = Yq(:,1:q-1)*pinv(Xq(:,1:q-1)); + obj.M = inv(Xq(:,1:q-1)*Xq(:,1:q-1)'); + end + obj.timestep = obj.timestep + q; + end + + function update(obj, xold, yold, xnew, ynew) + + % Update the DMD computation by sliding the finite time window forward + % Forget the oldest pair of snapshots (xold, yold), and includes 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 + % Xw = [x(k-w+2),x(k-w+3),...,x(k+1)], Yw = [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 + % xold = x(k-w+2), yold = y(k-w+2), xnew = x(k+1), ynew = y(k+1) + % Usage: wdmd.update(xold, yold, xnew, ynew) + + % Compute gamma + gamma = 1/(1+xnew'*(obj.M*xnew)); + % Compute Pk+1 + Pk1 = obj.M - gamma*((obj.M*xnew)*(xnew'*obj.M)); + % Compute beta + beta = 1/(1-xold'*(Pk1*xold)); + + % Update A + obj.A = obj.B + gamma*((ynew-obj.B*xnew)*(xnew'*obj.M)); + % Update B + obj.B = obj.A + beta*((-yold+obj.A*xold)*(xold'*Pk1)); + % Update M + obj.M = Pk1 + beta*((Pk1*xold)*(xold'*Pk1)); + + obj.timestep = obj.timestep + 1; + end + + function [evals, modes] = computemodes(obj) + % Compute and return DMD eigenvalues and DMD modes at current time step + % Usage: [evals, modes] = wdmd.computemodes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file diff --git a/matlab/online_demo.m b/matlab/online_demo.m new file mode 100644 index 0000000..50f45be --- /dev/null +++ b/matlab/online_demo.m @@ -0,0 +1,117 @@ +% An example to demonstrate online dynamic mode decomposition +% +% We take a 2D time varying system given by dx/dt = A(t)x +% where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +% w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +% are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit +% +% At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], +% that contain all the past snapshot pairs, we would like to compute +% Ak = Yk*pinv(Xk). This can be done by brute-force batch DMD, +% and by efficient rank-1 updating online DMD algrithm. +% +% Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly +% +% Online DMD computes the DMD matrix by using efficient rank-1 update idea +% +% We compare the performance of online DMD (with lambda=1,0.9) with the brute-force batch DMD +% approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Date created: April 2017 + +% define dynamics +epsilon = 1e-1; +dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x; +% generate data +tspan = [0 10]; +x0 = [1;0]; +[t,x] = ode45(dyn, tspan, x0); +% interpolate uniform time step +dt = 1e-1; +time = 0:dt:max(tspan); +xq = interp1(t,x,time); xq = xq'; +% extract snapshot pairs +x = xq(:,1:end-1); y = xq(:,2:end); t = time(2:end); +% true dynamics, eigenvalues +[n, m] = size(x); +A = zeros(n,n,m); +evals = zeros(n,m); +for k = 1:m + A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics + evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues +end + + +% visualize snapshots +figure, hold on +plot(time,xq(1,:),'x-',time,xq(2,:),'o-','LineWidth',2) +xlabel('Time','Interpreter','latex') +title('Snapshots','Interpreter','latex') +fl = legend('$x_1(t)$','$x_2(t)$'); +set(fl,'Interpreter','latex'); +box on +set(gca,'FontSize',20,'LineWidth',2) + + +% batch DMD +q = 20; +AbatchDMD = zeros(n,n,m); +evalsbatchDMD = zeros(n,m); +tic +for k = q+1:m + AbatchDMD(:,:,k) = y(:,1:k)*pinv(x(:,1:k)); + evalsbatchDMD(:,k) = log(eig(AbatchDMD(:,:,k)))/dt; +end +elapsed_time = toc; +fprintf('Batch DMD, elapsed time: %f seconds\n', elapsed_time) + +% Online DMD lambda = 1 +q = 20; +evalsonlineDMD1 = zeros(n,m); +% creat object and initialize with first q snapshot pairs +odmd = OnlineDMD(n,1); +odmd.initialize(x(:,1:q),y(:,1:q)); +% online DMD +tic +for k = q+1:m + odmd.update(x(:,k),y(:,k)); + evalsonlineDMD1(:,k) = log(eig(odmd.A))/dt; +end +elapsed_time = toc; +fprintf('Online DMD, lambda = 1, elapsed time: %f seconds\n', elapsed_time) + +% Online DMD, lambda = 0.9 +q = 20; +evalsonlineDMD09 = zeros(n,m); +% creat object and initialize with first q snapshot pairs +odmd = OnlineDMD(n,0.9); +odmd.initialize(x(:,1:q),y(:,1:q)); +% online DMD +tic +for k = q+1:m + odmd.update(x(:,k),y(:,k)); + evalsonlineDMD09(:,k) = log(eig(odmd.A))/dt; +end +elapsed_time = toc; +fprintf('Online DMD, lambda = 0.9, elapsed time: %f seconds\n', elapsed_time) + + +% visualize imaginary part of the continous time eigenvalues +% from true, batch, online (lambda=1), and online (lambda=0.9) +updateindex = q+1:m; +figure, hold on +plot(t,imag(evals(1,:)),'k-','LineWidth',3) +plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',3) +plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',3) +plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',3) +xlabel('Time','Interpreter','latex'), ylabel('Im') +title('Imaginary part of eigenvalues','Interpreter','latex') +fl = legend('True','batch','online, $\lambda=1$','online, $\lambda=0.9$'); +set(fl,'Interpreter','latex','Location','northwest'); +ylim([1,2]), xlim([0,10]) +box on +set(gca,'FontSize',18,'LineWidth',2) \ No newline at end of file diff --git a/matlab/window_demo.m b/matlab/window_demo.m new file mode 100644 index 0000000..e77770a --- /dev/null +++ b/matlab/window_demo.m @@ -0,0 +1,103 @@ +% An example to demonstrate window dynamic mode decomposition +% +% We take a 2D time varying system given by dx/dt = A(t)x +% where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +% w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +% are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit +% +% At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], +% that contain the recent w snapshot pairs from a finite time window, +% we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, +% and by efficient rank-2 updating window DMD algrithm. +% +% Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly +% +% Window DMD computes the DMD matrix by using efficient rank-2 update idea +% +% We compare the performance of window DMD with the brute-force mini-batch DMD +% approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% Date created: April 2017 + +% define dynamics +epsilon = 1e-1; +dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x; +% generate data +tspan = [0 10]; +x0 = [1;0]; +[t,x] = ode45(dyn, tspan, x0); +% interpolate uniform time step +dt = 1e-1; +time = 0:dt:max(tspan); +xq = interp1(t,x,time); xq = xq'; +% extract snapshot pairs +x = xq(:,1:end-1); y = xq(:,2:end); t = time(2:end); +% true dynamics, eigenvalues +[n, m] = size(x); +A = zeros(n,n,m); +evals = zeros(n,m); +for k = 1:m + A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics + evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues +end + + +% visualize snapshots +figure, hold on +plot(time,xq(1,:),'x-',time,xq(2,:),'o-','LineWidth',2) +xlabel('Time','Interpreter','latex') +title('Snapshots','Interpreter','latex') +fl = legend('$x_1(t)$','$x_2(t)$'); +set(fl,'Interpreter','latex'); +box on +set(gca,'FontSize',20,'LineWidth',2) + + +% mini-batch DMD +w = 20; % storage time window size, store recent w snapshot pairs +AminibatchDMD = zeros(n,n,m); +evalsminibatchDMD = zeros(n,m); +% mini-batch DMD +tic +for k = w+1:m + AminibatchDMD(:,:,k) = y(:,k-w+1:k)*pinv(x(:,k-w+1:k)); + evalsminibatchDMD(:,k) = log(eig(AminibatchDMD(:,:,k)))/dt; +end +elapsed_time = toc; +fprintf('Mini-batch DMD, elapsed time: %f seconds\n', elapsed_time) + + +% window DMD +w = 20; +evalswindowDMD = zeros(n,m); +% creat object and initialize with first w snapshot pairs +wdmd = WindowDMD(n,w); +wdmd.initialize(x(:,1:w), y(:,1:w)); +% window DMD +tic +for k = w+1:m + wdmd.update(x(:,k-w+1), y(:,k-w+1), x(:,k), y(:,k)); + evalswindowDMD(:,k) = log(eig(wdmd.A))/dt; +end +elapsed_time = toc; +fprintf('Window DMD, elapsed time: %f seconds\n', elapsed_time) + + +% visualize imaginary part of the continous time eigenvalues +% from true, mini-batch, and window +updateindex = w+1:m; +figure, hold on +plot(t,imag(evals(1,:)),'k-','LineWidth',3) +plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',3) +plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',3) +xlabel('Time','Interpreter','latex'), ylabel('Im') +title('Imaginary part of eigenvalues','Interpreter','latex') +fl = legend('True','mini-batch','window'); +set(fl,'Interpreter','latex','Location','northwest'); +ylim([1,2]), xlim([0,10]) +box on +set(gca,'FontSize',18,'LineWidth',2) \ No newline at end of file diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py new file mode 100644 index 0000000..28124fb --- /dev/null +++ b/python/dmdtools/online.py @@ -0,0 +1,109 @@ +import numpy as np + + +class OnlineDMD: + """OnlineDMD is a class that implements online dynamic mode decomposition + The time complexity (for one iteration) is O(n^2), and space complexity is + O(n^2), where n is the state dimension. + + Algorithm description: + At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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 correponding to consecutive states z(k-1) and z(k). + We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively + by efficient rank-1 updating online DMD algrithm. + + Usage: + odmd = OnlineDMD(n,forgetting) + odmd.initialize(Xq,Yq) + odmd.initilizeghost() + odmd.update(x,y) + evals, modes = odmd.computemodes() + + properties: + n: state dimension + forgetting: forgetting factor between (0,1] + timestep: number of snapshot pairs processed (i.e., the current time step) + 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 first 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 DMD computation when new snapshot pair (x,y) becomes available + computemodes(), compute and return DMD eigenvalues and DMD modes + + Authors: + Hao Zhang + Clarence W. Rowley + + Date created: April 2017 + + To import the OnlineDMD class, add import online at head of Python scripts. + To look up this documentation, type help(online.OnlineDMD) or online.OnlineDMD? + """ + def __init__(self, n=0, forgetting=1, timestep=0, A=None, P=None): + """ + Creat an object for online DMD + Usage: odmd = OnlineDMD(n,forgetting) + """ + self.n = n + self.forgetting = forgetting + self.timestep = timestep + if A is None or P is None: + self.A = np.zeros([n,n]) + self.P = np.zeros([n,n]) + else: + self.A = A + self.P = P + + def initialize(self, Xq, Yq): + """Initialize online DMD with first q snapshot pairs stored in (Xq, Yq) + Usage: odmd.initialize(Xq,Yq) + """ + q = len(Xq[0,:]) + Xqhat, Yqhat = np.zeros(Xq.shape), np.zeros(Yq.shape) + if self.timestep == 0 and self.n <= q: + sqrtlambda = np.sqrt(self.forgetting) + # multiply forgetting factor with snapshots + for i in range(q): + Xqhat[:,i] = Xq[:,i]*sqrtlambda**(q-1-i) + Yqhat[:,i] = Yq[:,i]*sqrtlambda**(q-1-i) + self.A = Yqhat.dot(np.linalg.pinv(Xqhat)) + self.P = np.linalg.inv(Xqhat.dot(Xqhat.T))/self.forgetting + self.timestep += q + + def initializeghost(self): + """Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs before t=0 + Usage: odmd.initilizeghost() + """ + epsilon=1e-15 + alpha = 1.0/epsilon + self.A = np.random.randn(self.n, self.n) + self.P = alpha*np.identity(self.n) + + def update(self, 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 = self.P.dot(x) + # compute gamma + gamma = 1.0/(1 + x.T.dot(Px)) + # update A + self.A += gamma*np.outer(y-self.A.dot(x),Px) + # update P + self.P = (self.P - gamma*np.outer(Px,Px))/self.forgetting + # time step + 1 + self.timestep += 1 + + def computemodes(self): + """Compute and return DMD eigenvalues and DMD modes at current time step + Usage: evals, modes = odmd.computemodes() + """ + evals, modes = np.linalg.eig(self.A) + return evals, modes \ No newline at end of file diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py new file mode 100644 index 0000000..bdaeeac --- /dev/null +++ b/python/dmdtools/window.py @@ -0,0 +1,112 @@ +import numpy as np + + +class WindowDMD: + """WindowDMD is a class that implements window dynamic mode decomposition + The time complexity (for one iteration) is O(n^2), and space complexity is + O(n^2), where n is the state dimension + + Algorithm description: + At time step k, define two matrix + Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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 correponding to consecutive states z(k-1) and z(k). + We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively + by efficient rank-2 updating window DMD algrithm. + + Usage: + wdmd = WindowDMD(n,windowsize) + wdmd.initialize(Xq,Yq) + wdmd.update(xold,yold,xnew,ynew) + evals, modes = wdmd.computemodes() + + properties: + n: state dimension + windowsize: window size + timestep: number of snapshot pairs processed (i.e., the current time step) + A: DMD matrix, size n by n + B: Intermediate DMD matrix for w-1 snapshot pairs, size n by n + M: Matrix that contains information about recent w-1 snapshots, size n by n + + methods: + initialize(Xq, Yq), initialize window DMD algorithm + update(x,y), update DMD computation when new snapshot pair (x,y) becomes available + computemodes(), compute and return DMD eigenvalues and DMD modes + + Authors: + Hao Zhang + Clarence W. Rowley + + Date created: April 2017 + + To import the WindowDMD class, add import window at head of Python scripts. + To look up this documentation, type help(window.WindowDMD) or window.WindowDMD? + """ + def __init__(self, n=0, windowsize=0, timestep=0, A=None, B=None, M=None): + """ + Creat an object for window DMD + Usage: wdmd = WindowDMD(n,windowsize) + """ + self.n = n + self.windowsize = windowsize + self.timestep = timestep + if A is None or B is None or M is None: + self.A = np.zeros([n,n]) + self.B = np.zeros([n,n]) + self.M = np.zeros([n,n]) + else: + self.A = A + self.B = B + self.M = M + + def initialize(self, Xq, Yq): + """Initialize window DMD with first q snapshot pairs stored in (Xq, Yq) + Usage: wdmd.initialize(Xq,Yq) + """ + q = len(Xq[0,:]) + if self.timestep == 0 and self.windowsize == q and self.windowsize >= self.n + 1: + self.A = Yq.dot(np.linalg.pinv(Xq)) + self.B = Yq[:,:q-1].dot(np.linalg.pinv(Xq[:,:q-1])) + self.M = np.linalg.inv(Xq[:,:q-1].dot(Xq[:,:q-1].T)) + self.timestep += q + + def update(self, xold, yold, xnew, ynew): + """Update the DMD computation by sliding the finite time window forward + Forget the oldest pair of snapshots (xold, yold), and includes 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 + xold = x(k-w+2), yold = y(k-w+2), xnew = x(k+1), ynew = y(k+1) + Usage: wdmd.update(xold, yold, xnew, ynew) + """ + # compute gamma + # compute M*xnew matrix vector product beforehand + Mxnew = self.M.dot(xnew) + gamma = 1.0/(1+xnew.T.dot(Mxnew)) + # compute Pk+1 + Pk1 = self.M - gamma*np.outer(Mxnew, Mxnew) + # compute beta + # compute P(k+1)*xold matrix vector product beforehand + Pk1xold = Pk1.dot(xold) + beta = 1.0/(1-xold.T.dot(Pk1xold)) + + # update A + self.A = self.B + gamma*np.outer(ynew - self.B.dot(xnew), Mxnew) + # update B + self.B = self.A + beta*np.outer(-yold + self.A.dot(xold), Pk1xold) + # update M + self.M = Pk1 + beta*np.outer(Pk1xold, Pk1xold) + + # time step + 1 + self.timestep += 1 + + def computemodes(self): + """Compute and return DMD eigenvalues and DMD modes at current time step + Usage: evals, modes = wdmd.computemodes() + """ + evals, modes = np.linalg.eig(self.A) + return evals, modes \ No newline at end of file diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py new file mode 100644 index 0000000..a1c261f --- /dev/null +++ b/python/scripts/online_demo.py @@ -0,0 +1,128 @@ +""" +An example to demonstrate online dynamic mode decomposition + +We take a 2D time varying system given by dx/dt = A(t)x +where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit + +At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], +that contain all the past snapshot pairs, we would like to compute +Ak = Yk*pinv(Xk). This can be done by brute-force batch DMD, +and by efficient rank-1 updating online DMD algrithm. + +Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly + +Online DMD computes the DMD matrix by using efficient rank-1 update idea + +We compare the performance of online DMD (with lambda=1,0.9) with the brute-force batch DMD +approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution + +Authors: + Hao Zhang + Clarence W. Rowley + +Date created: April 2017 +""" + + +import sys +sys.path.append('..') + +from online import OnlineDMD +import numpy as np +from scipy.integrate import odeint +import time +import matplotlib.pyplot as plt + + +# define dynamics +epsilon = 1e-1 +def dyn(x,t): + x1, x2 = x + dxdt = [(1+epsilon*t)*x2,-(1+epsilon*t)*x1] + return dxdt +# integrate from initial condition [1,0] +tspan = np.linspace(0,10,101) +dt = 0.1 +x0 = [1,0] +xsol = odeint(dyn,x0,tspan).T +x, y = xsol[:,:-1], xsol[:,1:] +t = tspan[:-1] +# true dynamics, true eigenvalues +n, m = len(x[:,0]), len(x[0,:]) +A = np.empty((n,n,m)) +evals = np.empty((n,m),dtype=complex) +for k in range(m): + A[:,:,k] = np.array([[0,(1+epsilon*t[k])],[-(1+epsilon*t[k]),0]]) + evals[:,k] = np.linalg.eigvals(A[:,:,k]) + + +# visualize snapshots +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(tspan, xsol[0,:], 'bs-', linewidth=2.0, label='$x_1(t)$') +plt.plot(tspan, xsol[1,:], 'g^-', linewidth=2.0, label='$x_2(t)$') +plt.legend(loc='best',fontsize=20 ,shadow=True) +plt.xlabel('Time', fontsize=20) +plt.title('Snapshots', fontsize=20) +plt.tick_params(labelsize=20) +plt.grid() +plt.show() + + +# batch DMD +q = 20 +AbatchDMD = np.empty((n,n,m)) +evalsbatchDMD = np.empty((n,m),dtype=complex) +start = time.clock() +for k in range(q,m): + AbatchDMD[:,:,k] = y[:,:k].dot(np.linalg.pinv(x[:,:k])) + evalsbatchDMD[:,k] = np.log(np.linalg.eigvals(AbatchDMD[:,:,k]))/dt +end = time.clock() +print "Batch DMD, time = " + str(end-start) + " secs" + + +# Online DMD, forgetting = 1 +q = 20 +evalsonlineDMD1 = np.empty((n,m),dtype=complex) +odmd = OnlineDMD(n,1.0) +odmd.initialize(x[:,:q],y[:,:q]) +start = time.clock() +for k in range(q,m): + odmd.update(x[:,k],y[:,k]) + evalsonlineDMD1[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt +end = time.clock() +print "Online DMD, forgetting = 1, time = " + str(end-start) + " secs" + + +# Online DMD, forgetting = 0.9 +q = 20 +evalsonlineDMD09 = np.empty((n,m),dtype=complex) +odmd = OnlineDMD(n,0.9) +odmd.initialize(x[:,:q],y[:,:q]) +start = time.clock() +for k in range(q,m): + odmd.update(x[:,k],y[:,k]) + evalsonlineDMD09[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt +end = time.clock() +print "Online DMD, forgetting = 0.9, time = " + str(end-start) + " secs" + + +# visualize true, batch, online (forgettting=1,0.9) +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='batch',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, $\lambda$=1',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, $\lambda$=0.9',linewidth=2.0) +plt.tick_params(labelsize=20) +plt.xlabel('Time', fontsize=20) +plt.ylabel('Im', fontsize=20) +plt.title('Imignary part of eigenvalues', fontsize=20) +plt.legend(loc='best', fontsize=20, shadow=True) +plt.xlim([0,10]) +plt.ylim([1,2]) +plt.show() \ No newline at end of file diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py new file mode 100644 index 0000000..de18b87 --- /dev/null +++ b/python/scripts/window_demo.py @@ -0,0 +1,114 @@ +""" +An example to demonstrate window dynamic mode decomposition + +We take a 2D time varying system given by dx/dt = A(t)x +where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit + +At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], +that contain the recent w snapshot pairs from a finite time window, +we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, +and by efficient rank-2 updating window DMD algrithm. + +Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly + +Window DMD computes the DMD matrix by using efficient rank-2 update idea + +We compare the performance of window DMD with the brute-force mini-batch DMD +approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution + +Authors: + Hao Zhang + Clarence W. Rowley + +Date created: April 2017 +""" + + +import sys +sys.path.append('..') + +from window import WindowDMD +import numpy as np +from scipy.integrate import odeint +import time +import matplotlib.pyplot as plt + + +# define dynamics +epsilon = 1e-1 +def dyn(x,t): + x1, x2 = x + dxdt = [(1+epsilon*t)*x2,-(1+epsilon*t)*x1] + return dxdt +# integrate from initial condition [1,0] +tspan = np.linspace(0,10,101) +dt = 0.1 +x0 = [1,0] +xsol = odeint(dyn,x0,tspan).T +x, y = xsol[:,:-1], xsol[:,1:] +t = tspan[:-1] +# true dynamics, true eigenvalues +n, m = len(x[:,0]), len(x[0,:]) +A = np.empty((n,n,m)) +evals = np.empty((n,m),dtype=complex) +for k in range(m): + A[:,:,k] = np.array([[0,(1+epsilon*t[k])],[-(1+epsilon*t[k]),0]]) + evals[:,k] = np.linalg.eigvals(A[:,:,k]) + + +# visualize snapshots +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(tspan, xsol[0,:], 'bs-', linewidth=2.0, label='$x_1(t)$') +plt.plot(tspan, xsol[1,:], 'g^-', linewidth=2.0, label='$x_2(t)$') +plt.legend(loc='best',fontsize=20 ,shadow=True) +plt.xlabel('Time', fontsize=20) +plt.title('Snapshots', fontsize=20) +plt.tick_params(labelsize=20) +plt.grid() +plt.show() + + +# mini-batch DMD, w = 20 +w = 20 +AminibatchDMD = np.empty((n,n,m)) +evalsminibatchDMD = np.empty((n,m),dtype=complex) +start = time.clock() +for k in range(w,m): + AminibatchDMD[:,:,k] = y[:,k-w:k].dot(np.linalg.pinv(x[:,k-w:k])) + evalsminibatchDMD[:,k] = np.log(np.linalg.eigvals(AminibatchDMD[:,:,k]))/dt +end = time.clock() +print "Mini-batch DMD, time = " + str(end-start) + " secs" + + +# Window DMD, w = 20 +w = 20 +evalswindowDMD = np.empty((n,m),dtype=complex) +wdmd = WindowDMD(n,w) +wdmd.initialize(x[:,:w],y[:,:w]) +start = time.clock() +for k in range(w,m): + wdmd.update(x[:,k-w+1],y[:,k-w+1],x[:,k],y[:,k]) + evalswindowDMD[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt +end = time.clock() +print "Window DMD, forgetting = 1, time = " + str(end-start) + " secs" + + +# visualize true, batch, window (forgettting=1,0.9) +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) +plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='mini-batch, w=20',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='window, w=20',linewidth=2.0) +plt.tick_params(labelsize=20) +plt.xlabel('Time', fontsize=20) +plt.ylabel('Im', fontsize=20) +plt.title('Imignary part of eigenvalues', fontsize=20) +plt.legend(loc='best', fontsize=20, shadow=True) +plt.xlim([0,10]) +plt.ylim([1,2]) +plt.show() \ No newline at end of file From 3018700bbe4e6371f90979c8e0d3e92afbd1e875 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 17:24:06 -0400 Subject: [PATCH 02/32] add readme for python online/window --- python/README.md | 42 +++++++++++++++++++++++++++++++++++ python/dmdtools/__init__.py | 4 ++-- python/dmdtools/online.py | 5 ++++- python/dmdtools/window.py | 5 ++++- python/scripts/online_demo.py | 6 ++--- python/scripts/window_demo.py | 4 ++-- 6 files changed, 57 insertions(+), 9 deletions(-) create mode 100644 python/README.md diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..cabcebf --- /dev/null +++ b/python/README.md @@ -0,0 +1,42 @@ +# README for a Online DMD and Window DMD +Python implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) + +## Online DMD algorithm description +At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. + +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Window DMD algorithm description +At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). + +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +by efficient rank-2 updating window DMD algroithm. +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Implementation +1.**online.py** implements **OnlineDMD** class in Python. +2.**window.py** implements **WindowDMD** class in Python. + +## Demos +1.**online_demo.py** demos the use of Python **OnlineDMD** class. +2.**window_demo.py** demos the use of Python **WindowDMD** class. + +## Authors: +Hao Zhang +Clarence W. Rowley + +## References +To be added + +## Date created: +April 2017 + diff --git a/python/dmdtools/__init__.py b/python/dmdtools/__init__.py index b6635bc..4aac6ee 100644 --- a/python/dmdtools/__init__.py +++ b/python/dmdtools/__init__.py @@ -1,4 +1,4 @@ from batch import * -#from regularized import * from streaming import * - +from online import * +from window import * \ No newline at end of file diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 28124fb..201ee50 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -35,9 +35,12 @@ class OnlineDMD: update(x,y), update DMD computation when new snapshot pair (x,y) becomes available computemodes(), compute and return DMD eigenvalues and DMD modes - Authors: + Authors: Hao Zhang Clarence W. Rowley + + References: + To be added. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index bdaeeac..38af65d 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -36,9 +36,12 @@ class WindowDMD: update(x,y), update DMD computation when new snapshot pair (x,y) becomes available computemodes(), compute and return DMD eigenvalues and DMD modes - Authors: + Authors: Hao Zhang Clarence W. Rowley + + References: + To be added. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index a1c261f..14bea41 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -29,7 +29,7 @@ import sys sys.path.append('..') -from online import OnlineDMD +import dmdtools import numpy as np from scipy.integrate import odeint import time @@ -87,7 +87,7 @@ def dyn(x,t): # Online DMD, forgetting = 1 q = 20 evalsonlineDMD1 = np.empty((n,m),dtype=complex) -odmd = OnlineDMD(n,1.0) +odmd = dmdtools.OnlineDMD(n,1.0) odmd.initialize(x[:,:q],y[:,:q]) start = time.clock() for k in range(q,m): @@ -100,7 +100,7 @@ def dyn(x,t): # Online DMD, forgetting = 0.9 q = 20 evalsonlineDMD09 = np.empty((n,m),dtype=complex) -odmd = OnlineDMD(n,0.9) +odmd = dmdtools.OnlineDMD(n,0.9) odmd.initialize(x[:,:q],y[:,:q]) start = time.clock() for k in range(q,m): diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index de18b87..3b4cb13 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -29,7 +29,7 @@ import sys sys.path.append('..') -from window import WindowDMD +import dmdtools import numpy as np from scipy.integrate import odeint import time @@ -87,7 +87,7 @@ def dyn(x,t): # Window DMD, w = 20 w = 20 evalswindowDMD = np.empty((n,m),dtype=complex) -wdmd = WindowDMD(n,w) +wdmd = dmdtools.WindowDMD(n,w) wdmd.initialize(x[:,:w],y[:,:w]) start = time.clock() for k in range(w,m): From 9882f930442a78f4b69124f20fa341c051e405f1 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 17:28:02 -0400 Subject: [PATCH 03/32] add readme for matlab online/window --- matlab/README.md | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 matlab/README.md diff --git a/matlab/README.md b/matlab/README.md new file mode 100644 index 0000000..d529f01 --- /dev/null +++ b/matlab/README.md @@ -0,0 +1,42 @@ +# README for a Online DMD and Window DMD +Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) + +## Online DMD algorithm description +At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. + +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Window DMD algorithm description +At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). + +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +by efficient rank-2 updating window DMD algroithm. +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Implementation +1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. +2.**WindomDMD.m** implements **WindowDMD** class in Matlab. + +## Demos +1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. +2.**window_demo.m** demos the use of Matlab **WindowDMD** class. + +## Authors: +Hao Zhang +Clarence W. Rowley + +## References +To be added + +## Date created: +April 2017 + From 80d70f6a8f2d2766bf308d3f810d5723ea100d01 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 17:35:18 -0400 Subject: [PATCH 04/32] update readme --- README.md | 2 +- python/README.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index e78e3fb..8401604 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,2 @@ # dmdtools -A library of tools for computing variants of Dynamic Mode Decomposition +A library of tools for computing variants of Dynamic Mode Decomposition \ No newline at end of file diff --git a/python/README.md b/python/README.md index cabcebf..712f02f 100644 --- a/python/README.md +++ b/python/README.md @@ -22,11 +22,11 @@ by efficient rank-2 updating window DMD algroithm. The time complexity (for one iteration) is O(n^2), and space complexity is O(n^2), where n is the state dimension. -## Implementation +## Implementation in folder **/dmdtools** 1.**online.py** implements **OnlineDMD** class in Python. 2.**window.py** implements **WindowDMD** class in Python. -## Demos +## Demos in folder **/scripts** 1.**online_demo.py** demos the use of Python **OnlineDMD** class. 2.**window_demo.py** demos the use of Python **WindowDMD** class. From dc512dd897f7e082d4c803fe813816620069eeaf Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 17:40:13 -0400 Subject: [PATCH 05/32] organize files --- matlab/{ => online}/OnlineDMD.m | 0 matlab/{ => online}/online_demo.m | 0 matlab/{ => streaming}/StreamingDMD.m | 0 matlab/{ => streaming}/StreamingTDMD.m | 0 matlab/{ => streaming}/sdmd_run.m | 0 matlab/{ => streaming}/stdmd_run.m | 0 matlab/{ => tdmd}/tdmd.m | 0 matlab/{ => tdmd}/tdmd_run.m | 0 matlab/{ => window}/WindowDMD.m | 0 matlab/{ => window}/window_demo.m | 0 10 files changed, 0 insertions(+), 0 deletions(-) rename matlab/{ => online}/OnlineDMD.m (100%) rename matlab/{ => online}/online_demo.m (100%) rename matlab/{ => streaming}/StreamingDMD.m (100%) rename matlab/{ => streaming}/StreamingTDMD.m (100%) rename matlab/{ => streaming}/sdmd_run.m (100%) rename matlab/{ => streaming}/stdmd_run.m (100%) rename matlab/{ => tdmd}/tdmd.m (100%) rename matlab/{ => tdmd}/tdmd_run.m (100%) rename matlab/{ => window}/WindowDMD.m (100%) rename matlab/{ => window}/window_demo.m (100%) diff --git a/matlab/OnlineDMD.m b/matlab/online/OnlineDMD.m similarity index 100% rename from matlab/OnlineDMD.m rename to matlab/online/OnlineDMD.m diff --git a/matlab/online_demo.m b/matlab/online/online_demo.m similarity index 100% rename from matlab/online_demo.m rename to matlab/online/online_demo.m diff --git a/matlab/StreamingDMD.m b/matlab/streaming/StreamingDMD.m similarity index 100% rename from matlab/StreamingDMD.m rename to matlab/streaming/StreamingDMD.m diff --git a/matlab/StreamingTDMD.m b/matlab/streaming/StreamingTDMD.m similarity index 100% rename from matlab/StreamingTDMD.m rename to matlab/streaming/StreamingTDMD.m diff --git a/matlab/sdmd_run.m b/matlab/streaming/sdmd_run.m similarity index 100% rename from matlab/sdmd_run.m rename to matlab/streaming/sdmd_run.m diff --git a/matlab/stdmd_run.m b/matlab/streaming/stdmd_run.m similarity index 100% rename from matlab/stdmd_run.m rename to matlab/streaming/stdmd_run.m diff --git a/matlab/tdmd.m b/matlab/tdmd/tdmd.m similarity index 100% rename from matlab/tdmd.m rename to matlab/tdmd/tdmd.m diff --git a/matlab/tdmd_run.m b/matlab/tdmd/tdmd_run.m similarity index 100% rename from matlab/tdmd_run.m rename to matlab/tdmd/tdmd_run.m diff --git a/matlab/WindowDMD.m b/matlab/window/WindowDMD.m similarity index 100% rename from matlab/WindowDMD.m rename to matlab/window/WindowDMD.m diff --git a/matlab/window_demo.m b/matlab/window/window_demo.m similarity index 100% rename from matlab/window_demo.m rename to matlab/window/window_demo.m From e53d2885c63249cd0517fac4d6586e1d0aab7800 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 17:52:50 -0400 Subject: [PATCH 06/32] organize files, update readme --- README.md | 15 ++++++++- matlab/README.md | 43 +++----------------------- matlab/online/README.md | 42 +++++++++++++++++++++++++ matlab/{README => streaming/README.md} | 0 matlab/window/README.md | 42 +++++++++++++++++++++++++ 5 files changed, 103 insertions(+), 39 deletions(-) create mode 100644 matlab/online/README.md rename matlab/{README => streaming/README.md} (100%) create mode 100644 matlab/window/README.md diff --git a/README.md b/README.md index 8401604..a0cf525 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,15 @@ # dmdtools -A library of tools for computing variants of Dynamic Mode Decomposition \ No newline at end of file +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 +4. Window DMD + +## algorithms implemented in Python +1. Standard batch processed DMD +2. Kernel DMD with polynomial kernel +3. Streaming DMD +4. Online DMD +5. Window DMD \ No newline at end of file diff --git a/matlab/README.md b/matlab/README.md index d529f01..95f84bd 100644 --- a/matlab/README.md +++ b/matlab/README.md @@ -1,42 +1,9 @@ -# README for a Online DMD and Window DMD -Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) +# README for Matlab implementation of variants of Dynamic Mode Decomposition -## Online DMD algorithm description -At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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. +## **tdmd** folder contains implementation of total-least-squares DMD (TDMD) and demo. -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). -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. +## **streaming** folder contains implementation of streaming (total-least-squares) DMD and demo. -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Window DMD algorithm description -At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). - -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively -by efficient rank-2 updating window DMD algroithm. -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Implementation -1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. -2.**WindomDMD.m** implements **WindowDMD** class in Matlab. - -## Demos -1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. -2.**window_demo.m** demos the use of Matlab **WindowDMD** class. - -## Authors: -Hao Zhang -Clarence W. Rowley - -## References -To be added - -## Date created: -April 2017 +## **online** folder contains implementation of online DMD and demo. +## **window** folder contains implementation of window DMD and demo. \ No newline at end of file diff --git a/matlab/online/README.md b/matlab/online/README.md new file mode 100644 index 0000000..d529f01 --- /dev/null +++ b/matlab/online/README.md @@ -0,0 +1,42 @@ +# README for a Online DMD and Window DMD +Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) + +## Online DMD algorithm description +At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. + +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Window DMD algorithm description +At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). + +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +by efficient rank-2 updating window DMD algroithm. +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Implementation +1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. +2.**WindomDMD.m** implements **WindowDMD** class in Matlab. + +## Demos +1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. +2.**window_demo.m** demos the use of Matlab **WindowDMD** class. + +## Authors: +Hao Zhang +Clarence W. Rowley + +## References +To be added + +## Date created: +April 2017 + diff --git a/matlab/README b/matlab/streaming/README.md similarity index 100% rename from matlab/README rename to matlab/streaming/README.md diff --git a/matlab/window/README.md b/matlab/window/README.md new file mode 100644 index 0000000..d529f01 --- /dev/null +++ b/matlab/window/README.md @@ -0,0 +1,42 @@ +# README for a Online DMD and Window DMD +Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) + +## Online DMD algorithm description +At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. + +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Window DMD algorithm description +At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). + +We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +by efficient rank-2 updating window DMD algroithm. +The time complexity (for one iteration) is O(n^2), and space complexity is +O(n^2), where n is the state dimension. + +## Implementation +1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. +2.**WindomDMD.m** implements **WindowDMD** class in Matlab. + +## Demos +1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. +2.**window_demo.m** demos the use of Matlab **WindowDMD** class. + +## Authors: +Hao Zhang +Clarence W. Rowley + +## References +To be added + +## Date created: +April 2017 + From ebfab70ca396b1b0f91ff2d8dfa28c7d8f24eff5 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 25 Apr 2017 18:01:15 -0400 Subject: [PATCH 07/32] update readme --- matlab/README.md | 12 +++++------ python/README.md | 54 +++++++++++------------------------------------- 2 files changed, 17 insertions(+), 49 deletions(-) diff --git a/matlab/README.md b/matlab/README.md index 95f84bd..178dd73 100644 --- a/matlab/README.md +++ b/matlab/README.md @@ -1,9 +1,7 @@ # README for Matlab implementation of variants of Dynamic Mode Decomposition -## **tdmd** folder contains implementation of total-least-squares DMD (TDMD) and demo. - -## **streaming** folder contains implementation of streaming (total-least-squares) DMD and demo. - -## **online** folder contains implementation of online DMD and demo. - -## **window** folder contains implementation of window DMD and demo. \ No newline at end of file +## implementations +1. **tdmd** folder contains implementation of total-least-squares DMD (TDMD) and demo. +2. **streaming** folder contains implementation of streaming DMD (including total-least-squares) and demo. +3. **online** folder contains implementation of online DMD and demo. +4. **window** folder contains implementation of window DMD and demo. \ No newline at end of file diff --git a/python/README.md b/python/README.md index 712f02f..0de879b 100644 --- a/python/README.md +++ b/python/README.md @@ -1,42 +1,12 @@ -# README for a Online DMD and Window DMD -Python implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) - -## Online DMD algorithm description -At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. - -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Window DMD algorithm description -At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). - -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively -by efficient rank-2 updating window DMD algroithm. -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Implementation in folder **/dmdtools** -1.**online.py** implements **OnlineDMD** class in Python. -2.**window.py** implements **WindowDMD** class in Python. - -## Demos in folder **/scripts** -1.**online_demo.py** demos the use of Python **OnlineDMD** class. -2.**window_demo.py** demos the use of Python **WindowDMD** class. - -## Authors: -Hao Zhang -Clarence W. Rowley - -## References -To be added - -## Date created: -April 2017 - +# README for Python implementation of variants of Dynamic Mode Decomposition + +## folders +1. **dmdtools** folder contains implementation of variants of Dynamic Mode Decomposition +2. **scripts** folder contains demos of various DMD algrithm +3. **tests** folder contains tests for regular DMD class and kernel DMD class. + +## implementations contained in **dmdtools** +1. **batch.py** implements Standard batch processed DMD and Kernel DMD with polynomial kernel +2. **streaming.py** implements Streaming DMD +3. **online.py** implements online DMD +4. **window.py** implements window DMD \ No newline at end of file From 78a9b7e7e9321beae84e7b8dc208b0fe8a3ba9fa Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 12 May 2017 12:32:20 -0400 Subject: [PATCH 08/32] keep original organization, update readme --- README.md | 4 +-- matlab/{online => }/OnlineDMD.m | 5 +++ matlab/README.md | 18 ++++++++--- matlab/{streaming => }/StreamingDMD.m | 0 matlab/{streaming => }/StreamingTDMD.m | 0 matlab/{window => }/WindowDMD.m | 5 +++ matlab/online/README.md | 42 -------------------------- matlab/{online => }/online_demo.m | 7 ++++- matlab/{streaming => }/sdmd_run.m | 0 matlab/{streaming => }/stdmd_run.m | 0 matlab/streaming/README.md | 32 -------------------- matlab/{tdmd => }/tdmd.m | 0 matlab/{tdmd => }/tdmd_run.m | 0 matlab/window/README.md | 42 -------------------------- matlab/{window => }/window_demo.m | 7 ++++- python/README.md | 18 +++++++---- python/dmdtools/online.py | 4 ++- python/dmdtools/window.py | 4 ++- python/scripts/online_demo.py | 7 ++++- python/scripts/window_demo.py | 5 +++ 20 files changed, 66 insertions(+), 134 deletions(-) rename matlab/{online => }/OnlineDMD.m (95%) rename matlab/{streaming => }/StreamingDMD.m (100%) rename matlab/{streaming => }/StreamingTDMD.m (100%) rename matlab/{window => }/WindowDMD.m (95%) delete mode 100644 matlab/online/README.md rename matlab/{online => }/online_demo.m (94%) rename matlab/{streaming => }/sdmd_run.m (100%) rename matlab/{streaming => }/stdmd_run.m (100%) delete mode 100755 matlab/streaming/README.md rename matlab/{tdmd => }/tdmd.m (100%) rename matlab/{tdmd => }/tdmd_run.m (100%) delete mode 100644 matlab/window/README.md rename matlab/{window => }/window_demo.m (93%) diff --git a/README.md b/README.md index a0cf525..1031a28 100644 --- a/README.md +++ b/README.md @@ -4,12 +4,12 @@ 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 +3. Online DMD (including weighted online DMD) 4. Window DMD ## algorithms implemented in Python 1. Standard batch processed DMD 2. Kernel DMD with polynomial kernel 3. Streaming DMD -4. Online DMD +4. Online DMD (including weighted online DMD) 5. Window DMD \ No newline at end of file diff --git a/matlab/online/OnlineDMD.m b/matlab/OnlineDMD.m similarity index 95% rename from matlab/online/OnlineDMD.m rename to matlab/OnlineDMD.m index a15772b..8603a30 100644 --- a/matlab/online/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -38,6 +38,11 @@ % Hao Zhang % Clarence W. Rowley % +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% in production, 2017. The paper will be available on arXiv soon. +% % Created: % April 2017. % diff --git a/matlab/README.md b/matlab/README.md index 178dd73..820dfea 100644 --- a/matlab/README.md +++ b/matlab/README.md @@ -1,7 +1,15 @@ # README for Matlab implementation of variants of Dynamic Mode Decomposition -## implementations -1. **tdmd** folder contains implementation of total-least-squares DMD (TDMD) and demo. -2. **streaming** folder contains implementation of streaming DMD (including total-least-squares) and demo. -3. **online** folder contains implementation of online DMD and demo. -4. **window** folder contains implementation of window DMD and demo. \ No newline at end of file +## 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. \ No newline at end of file diff --git a/matlab/streaming/StreamingDMD.m b/matlab/StreamingDMD.m similarity index 100% rename from matlab/streaming/StreamingDMD.m rename to matlab/StreamingDMD.m diff --git a/matlab/streaming/StreamingTDMD.m b/matlab/StreamingTDMD.m similarity index 100% rename from matlab/streaming/StreamingTDMD.m rename to matlab/StreamingTDMD.m diff --git a/matlab/window/WindowDMD.m b/matlab/WindowDMD.m similarity index 95% rename from matlab/window/WindowDMD.m rename to matlab/WindowDMD.m index b13c143..0456f65 100644 --- a/matlab/window/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -36,6 +36,11 @@ % Hao Zhang % Clarence W. Rowley % +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% in production, 2017. The paper will be available on arXiv soon. +% % Created: % April 2017. % diff --git a/matlab/online/README.md b/matlab/online/README.md deleted file mode 100644 index d529f01..0000000 --- a/matlab/online/README.md +++ /dev/null @@ -1,42 +0,0 @@ -# README for a Online DMD and Window DMD -Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) - -## Online DMD algorithm description -At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. - -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Window DMD algorithm description -At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). - -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively -by efficient rank-2 updating window DMD algroithm. -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Implementation -1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. -2.**WindomDMD.m** implements **WindowDMD** class in Matlab. - -## Demos -1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. -2.**window_demo.m** demos the use of Matlab **WindowDMD** class. - -## Authors: -Hao Zhang -Clarence W. Rowley - -## References -To be added - -## Date created: -April 2017 - diff --git a/matlab/online/online_demo.m b/matlab/online_demo.m similarity index 94% rename from matlab/online/online_demo.m rename to matlab/online_demo.m index 50f45be..ae721c6 100644 --- a/matlab/online/online_demo.m +++ b/matlab/online_demo.m @@ -20,7 +20,12 @@ % Authors: % Hao Zhang % Clarence W. Rowley -% +% +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% in production, 2017. The paper will be available on arXiv soon. +% % Date created: April 2017 % define dynamics diff --git a/matlab/streaming/sdmd_run.m b/matlab/sdmd_run.m similarity index 100% rename from matlab/streaming/sdmd_run.m rename to matlab/sdmd_run.m diff --git a/matlab/streaming/stdmd_run.m b/matlab/stdmd_run.m similarity index 100% rename from matlab/streaming/stdmd_run.m rename to matlab/stdmd_run.m diff --git a/matlab/streaming/README.md b/matlab/streaming/README.md deleted file mode 100755 index c448a34..0000000 --- a/matlab/streaming/README.md +++ /dev/null @@ -1,32 +0,0 @@ -README for a Streaming Dynamic Mode Decomposition (SDMD) -======================================================== - -The Matlab scripts contained in this distribution contain: (1) a Matlab class for an -object-oriented impelentation of incrementally updated DMD, and (2) an example -to demonstrate the use of this class on an arbitrary dynamical system. These scripts -are supplementary materials associated with Hemati et al. (2014). - - -SDMD Distribution Contents -========================== - -StreamingDMD.m --------------- -A Matlab class that implements SDMD in an object-oriented fashion, -based on the formulation in Hemati et al. (2014). - -sdmd_run.m ----------- -Demonstrates usage of the StreamingDMD class in solving an arbitrary -dynamical system with two characteristic frequencies. - - -Note: use the `help' and `doc' commands within Matlab to access more detailed -descriptions of these scripts. - - -References -========== -Maziar S. Hemati, Matthew O. Williams, and Clarence W. Rowley, -``Dynamic Mode Decomposition for Large and Streaming Datasets,'' -Physics of Fluids (2014). \ No newline at end of file diff --git a/matlab/tdmd/tdmd.m b/matlab/tdmd.m similarity index 100% rename from matlab/tdmd/tdmd.m rename to matlab/tdmd.m diff --git a/matlab/tdmd/tdmd_run.m b/matlab/tdmd_run.m similarity index 100% rename from matlab/tdmd/tdmd_run.m rename to matlab/tdmd_run.m diff --git a/matlab/window/README.md b/matlab/window/README.md deleted file mode 100644 index d529f01..0000000 --- a/matlab/window/README.md +++ /dev/null @@ -1,42 +0,0 @@ -# README for a Online DMD and Window DMD -Matlab implementation of online dynamic mode decomposition (Online DMD) and window dynamic mode decomposition (Window DMD) - -## Online DMD algorithm description -At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [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). -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algorithm. - -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Window DMD algorithm description -At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). - -We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively -by efficient rank-2 updating window DMD algroithm. -The time complexity (for one iteration) is O(n^2), and space complexity is -O(n^2), where n is the state dimension. - -## Implementation -1.**OnlineDMD.m** implements **OnlineDMD** class in Matlab. -2.**WindomDMD.m** implements **WindowDMD** class in Matlab. - -## Demos -1.**online_demo.m** demos the use of Matlab **OnlineDMD** class. -2.**window_demo.m** demos the use of Matlab **WindowDMD** class. - -## Authors: -Hao Zhang -Clarence W. Rowley - -## References -To be added - -## Date created: -April 2017 - diff --git a/matlab/window/window_demo.m b/matlab/window_demo.m similarity index 93% rename from matlab/window/window_demo.m rename to matlab/window_demo.m index e77770a..c4105dd 100644 --- a/matlab/window/window_demo.m +++ b/matlab/window_demo.m @@ -20,7 +20,12 @@ % Authors: % Hao Zhang % Clarence W. Rowley -% +% +% Reference: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% in production, 2017. The paper will be available on arXiv soon. +% % Date created: April 2017 % define dynamics diff --git a/python/README.md b/python/README.md index 0de879b..ac4f9a6 100644 --- a/python/README.md +++ b/python/README.md @@ -3,10 +3,16 @@ ## folders 1. **dmdtools** folder contains implementation of variants of Dynamic Mode Decomposition 2. **scripts** folder contains demos of various DMD algrithm -3. **tests** folder contains tests for regular DMD class and kernel DMD class. +3. **tests** folder contains tests for DMD class and kernel DMD class. -## implementations contained in **dmdtools** -1. **batch.py** implements Standard batch processed DMD and Kernel DMD with polynomial kernel -2. **streaming.py** implements Streaming DMD -3. **online.py** implements online DMD -4. **window.py** implements window DMD \ No newline at end of file +## implementations contained in **dmdtools** folder +1. **batch.py** implements Standard batch processed DMD (including total-least-squares) and Kernel DMD with polynomial kernel. +2. **streaming.py** implements Streaming DMD. +3. **online.py** implements online DMD (including weighted online DMD). +4. **window.py** implements window DMD. + +## demostrations contained in **scripts** folder +1. **total_dmd_example.py** demostrates the use of DMD and total-least-squares DMD. +2. **streaming_dmd_example.py** demostrates the use of Streaming DMD. +3. **online_demo.py** demostrates the use of online DMD +4. **window_demo.py** demostrates the use of window DMD \ No newline at end of file diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 201ee50..9797b95 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -40,7 +40,9 @@ class OnlineDMD: Clarence W. Rowley References: - To be added. + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + in production, 2017. The paper will be available on arXiv soon. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 38af65d..790ac3e 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -41,7 +41,9 @@ class WindowDMD: Clarence W. Rowley References: - To be added. + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + in production, 2017. The paper will be available on arXiv soon. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 14bea41..54ded77 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -21,7 +21,12 @@ Authors: Hao Zhang Clarence W. Rowley - + +References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + in production, 2017. The paper will be available on arXiv soon. + Date created: April 2017 """ diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 3b4cb13..83ffdcf 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -22,6 +22,11 @@ Hao Zhang Clarence W. Rowley +References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + in production, 2017. The paper will be available on arXiv soon. + Date created: April 2017 """ From 9db5a73e8be50948f87442c80e37a7559daf2a3b Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Sat, 13 May 2017 16:26:34 -0400 Subject: [PATCH 09/32] update reference info --- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 2 +- matlab/online_demo.m | 2 +- matlab/window_demo.m | 2 +- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 2 +- python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 8603a30..34264ab 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -41,7 +41,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", -% in production, 2017. The paper will be available on arXiv soon. +% in production, 2017. To be submitted for publication, available on arXiv. % % Created: % April 2017. diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 0456f65..ba97959 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -39,7 +39,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", -% in production, 2017. The paper will be available on arXiv soon. +% in production, 2017. To be submitted for publication, available on arXiv. % % Created: % April 2017. diff --git a/matlab/online_demo.m b/matlab/online_demo.m index ae721c6..57e51ed 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -24,7 +24,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", -% in production, 2017. The paper will be available on arXiv soon. +% in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/matlab/window_demo.m b/matlab/window_demo.m index c4105dd..053d200 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -24,7 +24,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", -% in production, 2017. The paper will be available on arXiv soon. +% in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 9797b95..e49689b 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -42,7 +42,7 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", - in production, 2017. The paper will be available on arXiv soon. + in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 790ac3e..734f566 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -43,7 +43,7 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", - in production, 2017. The paper will be available on arXiv soon. + in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 54ded77..dcf24e1 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -25,7 +25,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", - in production, 2017. The paper will be available on arXiv soon. + in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 """ diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 83ffdcf..1934333 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -25,7 +25,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", - in production, 2017. The paper will be available on arXiv soon. + in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 """ From a53e63fabd00c65918bd5aa19edf7320f604707e Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Sun, 14 May 2017 14:05:38 -0400 Subject: [PATCH 10/32] update reference info --- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 2 +- matlab/online_demo.m | 2 +- matlab/window_demo.m | 2 +- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 2 +- python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 34264ab..455ac17 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -40,7 +40,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", % in production, 2017. To be submitted for publication, available on arXiv. % % Created: diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index ba97959..d29555e 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -38,7 +38,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", % in production, 2017. To be submitted for publication, available on arXiv. % % Created: diff --git a/matlab/online_demo.m b/matlab/online_demo.m index 57e51ed..09a12bc 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -23,7 +23,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", % in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 053d200..3906065 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -23,7 +23,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", +% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", % in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index e49689b..757537e 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -41,7 +41,7 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 734f566..1e38490 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -42,7 +42,7 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index dcf24e1..25573ab 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -24,7 +24,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 1934333..2e31545 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -24,7 +24,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic Time Methods for Online Dynamic Mode Decomposition", + ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 From 6076eede0e439a1fd9c393e6bef1fee119acab4b Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Mon, 15 May 2017 10:52:24 -0400 Subject: [PATCH 11/32] update online DMD reference title --- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 2 +- matlab/online_demo.m | 2 +- matlab/window_demo.m | 2 +- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 2 +- python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 455ac17..66f410b 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -40,7 +40,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", +% ``Online Dynamic Mode Decomposition for Time-varying Systems", % in production, 2017. To be submitted for publication, available on arXiv. % % Created: diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index d29555e..ca47e43 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -38,7 +38,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", +% ``Online Dynamic Mode Decomposition for Time-varying Systems", % in production, 2017. To be submitted for publication, available on arXiv. % % Created: diff --git a/matlab/online_demo.m b/matlab/online_demo.m index 09a12bc..30c3802 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -23,7 +23,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", +% ``Online Dynamic Mode Decomposition for Time-varying Systems", % in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 3906065..4323a31 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -23,7 +23,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", +% ``Online Dynamic Mode Decomposition for Time-varying Systems", % in production, 2017. To be submitted for publication, available on arXiv. % % Date created: April 2017 diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 757537e..ef4a778 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -41,7 +41,7 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", + ``Online Dynamic Mode Decomposition for Time-varying Systems", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 1e38490..7f2ba01 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -42,7 +42,7 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", + ``Online Dynamic Mode Decomposition for Time-varying Systems", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 25573ab..b4041da 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -24,7 +24,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", + ``Online Dynamic Mode Decomposition for Time-varying Systems", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 2e31545..41d7ed3 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -24,7 +24,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Fast Quadratic-time Methods for Online Dynamic Mode Decomposition", + ``Online Dynamic Mode Decomposition for Time-varying Systems", in production, 2017. To be submitted for publication, available on arXiv. Date created: April 2017 From fbca7579c15505d5eedf795d1f1e2db4a09ae4c3 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Thu, 18 May 2017 13:36:14 -0400 Subject: [PATCH 12/32] make matlabe online/window more efficient --- matlab/OnlineDMD.m | 18 +++++++++--------- matlab/WindowDMD.m | 17 +++++++++++------ matlab/online_demo.m | 4 ++-- matlab/window_demo.m | 4 ++-- 4 files changed, 24 insertions(+), 19 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 66f410b..5f4bb5e 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -74,11 +74,9 @@ function initialize(obj, Xq, Yq) % Usage: odmd.initialize(Xq,Yq) q = length(Xq(1,:)); if(obj.timestep == 0 && q>=obj.n) - sqrtlambda = sqrt(obj.lambda); - for i = 1:q - Xq(:,i) = Xq(:,i)*sqrtlambda^(q-i); - Yq(:,i) = Yq(:,i)*sqrtlambda^(q-i); - end + weight = (sqrt(obj.lambda)).^(q-1:-1:0); + Xq = Xq.*weight; + Yq = Yq.*weight; obj.A = Yq*pinv(Xq); obj.P = inv(Xq*Xq')/obj.lambda; end @@ -100,13 +98,15 @@ function update(obj, 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'*(obj.P*x)); + gamma = 1/(1+x'*Px); % Update A - obj.A = obj.A + gamma*((y-obj.A*x)*(x'*obj.P)); + obj.A = obj.A + gamma*((y-obj.A*x)*Px'); % Update P - obj.P = (obj.P - gamma*((obj.P*x)*(x'*obj.P)))/obj.lambda; - + obj.P = (obj.P - gamma*(Px*Px'))/obj.lambda; + % time step + 1 obj.timestep = obj.timestep + 1; end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index ca47e43..a168e06 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -92,20 +92,25 @@ function update(obj, xold, yold, xnew, ynew) % xold = x(k-w+2), yold = y(k-w+2), xnew = x(k+1), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) + % compute M*xnew matrix vector product beforehand + Mxnew = obj.M*xnew; % Compute gamma - gamma = 1/(1+xnew'*(obj.M*xnew)); + gamma = 1/(1+xnew'*(Mxnew)); % Compute Pk+1 - Pk1 = obj.M - gamma*((obj.M*xnew)*(xnew'*obj.M)); + Pk1 = obj.M - gamma*(Mxnew*Mxnew'); + % compute P(k+1)*xold matrix vector product beforehand + Pk1xold = Pk1*xold; % Compute beta - beta = 1/(1-xold'*(Pk1*xold)); + beta = 1/(1-xold'*(Pk1xold)); % Update A - obj.A = obj.B + gamma*((ynew-obj.B*xnew)*(xnew'*obj.M)); + obj.A = obj.B + gamma*((ynew-obj.B*xnew)*Mxnew'); % Update B - obj.B = obj.A + beta*((-yold+obj.A*xold)*(xold'*Pk1)); + obj.B = obj.A + beta*((-yold+obj.A*xold)*Pk1xold'); % Update M - obj.M = Pk1 + beta*((Pk1*xold)*(xold'*Pk1)); + obj.M = Pk1 + beta*(Pk1xold*Pk1xold'); + % time step + 1 obj.timestep = obj.timestep + 1; end diff --git a/matlab/online_demo.m b/matlab/online_demo.m index 30c3802..0cfd4c7 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -113,10 +113,10 @@ plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',3) plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',3) plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',3) -xlabel('Time','Interpreter','latex'), ylabel('Im') +xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') title('Imaginary part of eigenvalues','Interpreter','latex') fl = legend('True','batch','online, $\lambda=1$','online, $\lambda=0.9$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on -set(gca,'FontSize',18,'LineWidth',2) \ No newline at end of file +set(gca,'FontSize',20,'LineWidth',2) \ No newline at end of file diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 4323a31..eb167a8 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -99,10 +99,10 @@ plot(t,imag(evals(1,:)),'k-','LineWidth',3) plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',3) plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',3) -xlabel('Time','Interpreter','latex'), ylabel('Im') +xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') title('Imaginary part of eigenvalues','Interpreter','latex') fl = legend('True','mini-batch','window'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on -set(gca,'FontSize',18,'LineWidth',2) \ No newline at end of file +set(gca,'FontSize',20,'LineWidth',2) \ No newline at end of file From b007fa471d0b7a62b9ba13a0cdfc9c879b554ddb Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Thu, 18 May 2017 22:24:31 -0400 Subject: [PATCH 13/32] update window implementation with A and P --- matlab/OnlineDMD.m | 4 +- matlab/WindowDMD.m | 59 ++++++++-------- matlab/WindowDMD.m~ | 128 ++++++++++++++++++++++++++++++++++ matlab/online_demo.m | 8 +-- matlab/window_demo.m | 8 +-- python/dmdtools/online.py | 4 +- python/dmdtools/window.py | 50 ++++++------- python/scripts/window_demo.py | 4 +- 8 files changed, 194 insertions(+), 71 deletions(-) create mode 100644 matlab/WindowDMD.m~ diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 5f4bb5e..8f03483 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -103,9 +103,9 @@ function update(obj, x, y) % Compute gamma gamma = 1/(1+x'*Px); % Update A - obj.A = obj.A + gamma*((y-obj.A*x)*Px'); + obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; % Update P - obj.P = (obj.P - gamma*(Px*Px'))/obj.lambda; + obj.P = (obj.P - (gamma*Px)*Px')/obj.lambda; % time step + 1 obj.timestep = obj.timestep + 1; end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index a168e06..c37a662 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -23,13 +23,15 @@ % n: state dimension % w: finite time window size % timestep: number of snapshot pairs processed -% A: Intermediate DMD matrix for w-1 snapshot pairs, size n by n -% B: Intermediate DMD matrix for w-1 snapshot pairs, size n by n -% M: Matrix that contains information about recent w-1 snapshots, size n by n +% 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(Xq, Yq), initialize window DMD algorithm -% update(xold, yold, xnew, ynew), update when new snapshot pair becomes available +% update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, +% and remeber new snapshot pair +% At time k+1, if 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 xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % computemodes(), compute and return DMD eigenvalues and DMD mdoes % % Authors: @@ -51,9 +53,8 @@ n = 0; % state dimension w = 0; % weighting factor timestep = 0; % number of snapshots processed - A; % Intermediate DMD matrix for w-1 snapshot pairs, size n by n - B; % Intermediate DMD matrix for w-1 snapshot pairs, size n by n - M; % Matrix that contains information about recent w-1 snapshots, size n by n + 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 @@ -64,8 +65,7 @@ obj.n = n; obj.w = w; obj.A = zeros(n,n); - obj.B = zeros(n,n); - obj.M = zeros(n,n); + obj.P = zeros(n,n); end end @@ -75,8 +75,7 @@ function initialize(obj, Xq, Yq) q = length(Xq(1,:)); if(obj.timestep == 0 && obj.w == q && obj.w >= obj.n+1) obj.A = Yq*pinv(Xq); - obj.B = Yq(:,1:q-1)*pinv(Xq(:,1:q-1)); - obj.M = inv(Xq(:,1:q-1)*Xq(:,1:q-1)'); + obj.P = inv(Xq*Xq'); end obj.timestep = obj.timestep + q; end @@ -84,31 +83,31 @@ function initialize(obj, Xq, Yq) function update(obj, xold, yold, xnew, ynew) % Update the DMD computation by sliding the finite time window forward - % Forget the oldest pair of snapshots (xold, yold), and includes the newest + % Forget the oldest pair of snapshots (xold, yold), and remebers 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 % Xw = [x(k-w+2),x(k-w+3),...,x(k+1)], Yw = [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 - % xold = x(k-w+2), yold = y(k-w+2), xnew = x(k+1), ynew = y(k+1) + % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) - % compute M*xnew matrix vector product beforehand - Mxnew = obj.M*xnew; - % Compute gamma - gamma = 1/(1+xnew'*(Mxnew)); - % Compute Pk+1 - Pk1 = obj.M - gamma*(Mxnew*Mxnew'); - % compute P(k+1)*xold matrix vector product beforehand - Pk1xold = Pk1*xold; - % Compute beta - beta = 1/(1-xold'*(Pk1xold)); - - % Update A - obj.A = obj.B + gamma*((ynew-obj.B*xnew)*Mxnew'); - % Update B - obj.B = obj.A + beta*((-yold+obj.A*xold)*Pk1xold'); - % Update M - obj.M = Pk1 + beta*(Pk1xold*Pk1xold'); + % compute Pk*xold matrix vector product beforehand + Pkxold = obj.P*xold; + % compute beta + beta = 1/(1-xold'*Pkxold); + % compute Ako + Ako = obj.A + (beta*(-yold+obj.A*xold))*Pkxold'; + % compute Pko + Pko = obj.P + (beta*Pkxold)*Pkxold'; + + % compute Pko*xnew matrix vector product beforehand + Pkoxnew = Pko*xnew; + % compute gamma + gamma = 1/(1+xnew'*Pkoxnew); + % update A + obj.A = Ako + (gamma*(ynew - Ako*xnew))*Pkoxnew'; + % update P + obj.P = Pko - (gamma*Pkoxnew)*Pkoxnew'; % time step + 1 obj.timestep = obj.timestep + 1; diff --git a/matlab/WindowDMD.m~ b/matlab/WindowDMD.m~ new file mode 100644 index 0000000..3161379 --- /dev/null +++ b/matlab/WindowDMD.m~ @@ -0,0 +1,128 @@ +% WindowDMD is a class that implements window dynamic mode decomposition +% The time complexity (for one iteration) is O(n^2), and space complexity is +% O(n^2), where n is the state dimension +% +% Algorithm description: +% At time step k, define two matrix +% Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-2 updating window DMD algorithm. +% +% Usage: +% wdmd = WindowDMD(n,w) +% wdmd.initialize(Xq,Yq) +% wdmd.update(xold, yold, xnew, ynew) +% [evals, modes] = wdmd.computemodes() +% +% properties: +% n: state dimension +% w: finite time window size +% timestep: number of snapshot pairs processed +% 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(Xq, Yq), initialize window DMD algorithm +% update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, +% and remeber new snapshot pair +% computemodes(), compute and return DMD eigenvalues and DMD mdoes +% +% 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", +% in production, 2017. To be submitted for publication, available on arXiv. +% +% Created: +% April 2017. +% +% To look up the documentation, type help WindowDMD + +classdef WindowDMD < handle + properties + n = 0; % state dimension + w = 0; % weighting factor + timestep = 0; % number of snapshots processed + 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) + % Creat an object for window DMD + % Usage: wdmd = WindowDMD(n,w) + if nargin == 2 + obj.n = n; + obj.w = w; + obj.A = zeros(n,n); + obj.P = zeros(n,n); + end + end + + function initialize(obj, Xq, Yq) + % Initialize WnlineDMD with q snapshot pairs stored in (Xq, Yq) + % Usage: wdmd.initialize(Xq,Yq) + q = length(Xq(1,:)); + if(obj.timestep == 0 && obj.w == q && obj.w >= obj.n+1) + obj.A = Yq*pinv(Xq); + obj.P = inv(Xq*Xq'); + end + obj.timestep = obj.timestep + q; + end + + function update(obj, xold, yold, xnew, ynew) + + % Update the DMD computation by sliding the finite time window forward + % Forget the oldest pair of snapshots (xold, yold), and remebers 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 + % Xw = [x(k-w+2),x(k-w+3),...,x(k+1)], Yw = [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 + % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) + % Usage: wdmd.update(xold, yold, xnew, ynew) + + % compute Pk*xold matrix vector product beforehand + Pxold = obj.P*xold; + % compute beta + beta = 1/(1-xold'*Pkxold); + % compute Pko + Pko = obj.P + (beta*Pkxold)*Pkxold'; + + + % compute M*xnew matrix vector product beforehand + Mxnew = obj.M*xnew; + % Compute gamma + gamma = 1/(1+xnew'*(Mxnew)); + % Compute Pk+1 + Pk1 = obj.M - gamma*(Mxnew*Mxnew'); + % compute P(k+1)*xold matrix vector product beforehand + Pk1xold = Pk1*xold; + % Compute beta + beta = 1/(1-xold'*(Pk1xold)); + + % Update A + obj.A = obj.B + gamma*((ynew-obj.B*xnew)*Mxnew'); + % Update B + obj.B = obj.A + beta*((-yold+obj.A*xold)*Pk1xold'); + % Update M + obj.M = Pk1 + beta*(Pk1xold*Pk1xold'); + + % time step + 1 + obj.timestep = obj.timestep + 1; + end + + function [evals, modes] = computemodes(obj) + % Compute and return DMD eigenvalues and DMD modes at current time step + % Usage: [evals, modes] = wdmd.computemodes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file diff --git a/matlab/online_demo.m b/matlab/online_demo.m index 0cfd4c7..b8f618d 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -109,10 +109,10 @@ % from true, batch, online (lambda=1), and online (lambda=0.9) updateindex = q+1:m; figure, hold on -plot(t,imag(evals(1,:)),'k-','LineWidth',3) -plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',3) -plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',3) -plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',3) +plot(t,imag(evals(1,:)),'k-','LineWidth',2) +plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) +plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) +plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') title('Imaginary part of eigenvalues','Interpreter','latex') fl = legend('True','batch','online, $\lambda=1$','online, $\lambda=0.9$'); diff --git a/matlab/window_demo.m b/matlab/window_demo.m index eb167a8..a824a2b 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -85,7 +85,7 @@ % window DMD tic for k = w+1:m - wdmd.update(x(:,k-w+1), y(:,k-w+1), x(:,k), y(:,k)); + wdmd.update(x(:,k-w), y(:,k-w), x(:,k), y(:,k)); evalswindowDMD(:,k) = log(eig(wdmd.A))/dt; end elapsed_time = toc; @@ -96,9 +96,9 @@ % from true, mini-batch, and window updateindex = w+1:m; figure, hold on -plot(t,imag(evals(1,:)),'k-','LineWidth',3) -plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',3) -plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',3) +plot(t,imag(evals(1,:)),'k-','LineWidth',2) +plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) +plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') title('Imaginary part of eigenvalues','Interpreter','latex') fl = legend('True','mini-batch','window'); diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index ef4a778..ae9ba8e 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -100,9 +100,9 @@ def update(self, x, y): # compute gamma gamma = 1.0/(1 + x.T.dot(Px)) # update A - self.A += gamma*np.outer(y-self.A.dot(x),Px) + self.A += np.outer(y-self.A.dot(x),gamma*Px) # update P - self.P = (self.P - gamma*np.outer(Px,Px))/self.forgetting + self.P = (self.P - np.outer(gamma*Px,Px))/self.forgetting # time step + 1 self.timestep += 1 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 7f2ba01..8d3350c 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -28,8 +28,7 @@ class WindowDMD: windowsize: window size timestep: number of snapshot pairs processed (i.e., the current time step) A: DMD matrix, size n by n - B: Intermediate DMD matrix for w-1 snapshot pairs, size n by n - M: Matrix that contains information about recent w-1 snapshots, size n by n + P: Matrix that contains information about recent w snapshots, size n by n methods: initialize(Xq, Yq), initialize window DMD algorithm @@ -50,7 +49,7 @@ class WindowDMD: To import the WindowDMD class, add import window at head of Python scripts. To look up this documentation, type help(window.WindowDMD) or window.WindowDMD? """ - def __init__(self, n=0, windowsize=0, timestep=0, A=None, B=None, M=None): + def __init__(self, n=0, windowsize=0, timestep=0, A=None, P=None): """ Creat an object for window DMD Usage: wdmd = WindowDMD(n,windowsize) @@ -58,14 +57,12 @@ def __init__(self, n=0, windowsize=0, timestep=0, A=None, B=None, M=None): self.n = n self.windowsize = windowsize self.timestep = timestep - if A is None or B is None or M is None: + if A is None or P is None: self.A = np.zeros([n,n]) - self.B = np.zeros([n,n]) - self.M = np.zeros([n,n]) + self.P = np.zeros([n,n]) else: self.A = A - self.B = B - self.M = M + self.P = P def initialize(self, Xq, Yq): """Initialize window DMD with first q snapshot pairs stored in (Xq, Yq) @@ -74,37 +71,36 @@ def initialize(self, Xq, Yq): q = len(Xq[0,:]) if self.timestep == 0 and self.windowsize == q and self.windowsize >= self.n + 1: self.A = Yq.dot(np.linalg.pinv(Xq)) - self.B = Yq[:,:q-1].dot(np.linalg.pinv(Xq[:,:q-1])) - self.M = np.linalg.inv(Xq[:,:q-1].dot(Xq[:,:q-1].T)) + self.P = np.linalg.inv(Xq.dot(Xq.T)) self.timestep += q def update(self, xold, yold, xnew, ynew): """Update the DMD computation by sliding the finite time window forward - Forget the oldest pair of snapshots (xold, yold), and includes the newest + 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 - xold = x(k-w+2), yold = y(k-w+2), xnew = x(k+1), ynew = y(k+1) + xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) Usage: wdmd.update(xold, yold, xnew, ynew) """ - # compute gamma - # compute M*xnew matrix vector product beforehand - Mxnew = self.M.dot(xnew) - gamma = 1.0/(1+xnew.T.dot(Mxnew)) - # compute Pk+1 - Pk1 = self.M - gamma*np.outer(Mxnew, Mxnew) + # compute Pkxold matrix vector product beforehand + Pkxold = self.P.dot(xold) # compute beta - # compute P(k+1)*xold matrix vector product beforehand - Pk1xold = Pk1.dot(xold) - beta = 1.0/(1-xold.T.dot(Pk1xold)) + beta = 1.0/(1-xold.T.dot(Pkxold)) + # compute Ako + Ako = self.A + np.outer(beta*(-yold+self.A.dot(xold)),Pkxold) + # compute Pko + Pko = self.P + np.outer(beta*Pkxold, Pkxold) - # update A - self.A = self.B + gamma*np.outer(ynew - self.B.dot(xnew), Mxnew) - # update B - self.B = self.A + beta*np.outer(-yold + self.A.dot(xold), Pk1xold) - # update M - self.M = Pk1 + beta*np.outer(Pk1xold, Pk1xold) + # compute Pko*xnew matrix vector product beforehand + Pkoxnew = Pko.dot(xnew) + # compute gamma + gamma = 1.0/(1+xnew.T.dot(Pkoxnew)) + # update Ak + self.A = Ako + np.outer(gamma*(ynew-Ako.dot(xnew)),Pkoxnew) + # update Pk + self.P = Pko - np.outer(gamma*Pkoxnew, Pkoxnew) # time step + 1 self.timestep += 1 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 41d7ed3..f78012a 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -96,10 +96,10 @@ def dyn(x,t): wdmd.initialize(x[:,:w],y[:,:w]) start = time.clock() for k in range(w,m): - wdmd.update(x[:,k-w+1],y[:,k-w+1],x[:,k],y[:,k]) + wdmd.update(x[:,k-w],y[:,k-w],x[:,k],y[:,k]) evalswindowDMD[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt end = time.clock() -print "Window DMD, forgetting = 1, time = " + str(end-start) + " secs" +print "Window DMD, time = " + str(end-start) + " secs" # visualize true, batch, window (forgettting=1,0.9) From 1b7e1d8751fd535e750a5867f21e361271ddf7aa Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 09:14:22 -0400 Subject: [PATCH 14/32] update weighting notation --- matlab/OnlineDMD.m | 25 ++++--- matlab/WindowDMD.m | 8 ++- matlab/WindowDMD.m~ | 128 ---------------------------------- matlab/online_demo.m | 19 ++--- matlab/window_demo.m | 5 +- python/dmdtools/online.py | 19 ++--- python/dmdtools/window.py | 4 ++ python/scripts/online_demo.py | 16 ++--- python/scripts/window_demo.py | 4 +- 9 files changed, 57 insertions(+), 171 deletions(-) delete mode 100644 matlab/WindowDMD.m~ diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 8f03483..dd998bb 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -3,16 +3,17 @@ % O(n^2), where n is the state dimension. % % Algorithm description: -% At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], +% 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. % % Usage: -% odmd = OnlineDMD(n,lambda) +% odmd = OnlineDMD(n,alpha) % odmd.initialize(Xq,Yq) % odmd.initilizeghost() % odmd.update(x,y) @@ -20,7 +21,7 @@ % % properties: % n: state dimension -% lambda: weighting factor between 0 and 1 +% alpha: weighting factor between 0 and 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 @@ -51,19 +52,19 @@ classdef OnlineDMD < handle properties n = 0; % state dimension - lambda = 1; % weighting factor + alpha = 1; % weighting factor timestep = 0; % number of snapshots processed A; % DMD matrix P; % matrix that contains information about past snapshots end methods - function obj = OnlineDMD(n,lambda) + function obj = OnlineDMD(n,alpha) % Creat an object for online DMD - % Usage: odmd = OnlineDMD(n,lambda) + % Usage: odmd = OnlineDMD(n,alpha) if nargin == 2 obj.n = n; - obj.lambda = lambda; + obj.alpha = alpha; obj.A = zeros(n,n); obj.P = zeros(n,n); end @@ -74,11 +75,11 @@ function initialize(obj, Xq, Yq) % Usage: odmd.initialize(Xq,Yq) q = length(Xq(1,:)); if(obj.timestep == 0 && q>=obj.n) - weight = (sqrt(obj.lambda)).^(q-1:-1:0); + weight = (sqrt(obj.alpha)).^(q-1:-1:0); Xq = Xq.*weight; Yq = Yq.*weight; obj.A = Yq*pinv(Xq); - obj.P = inv(Xq*Xq')/obj.lambda; + obj.P = inv(Xq*Xq')/obj.alpha; end obj.timestep = obj.timestep + q; end @@ -86,10 +87,8 @@ function initialize(obj, Xq, Yq) function initializeghost(obj) % Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs before t=0 % Usage: odmd.initilizeghost() - epsilon = 1e-15; - alpha = 1.0/epsilon; obj.A = randn(obj.n, obj.n); - obj.P = alpha*eye(obj.n); + obj.P = (1/eps)*eye(obj.n); end function update(obj, x, y) @@ -105,7 +104,7 @@ function update(obj, x, y) % Update A obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; % Update P - obj.P = (obj.P - (gamma*Px)*Px')/obj.lambda; + obj.P = (obj.P - (gamma*Px)*Px')/obj.alpha; % time step + 1 obj.timestep = obj.timestep + 1; end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index c37a662..e81713d 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -4,12 +4,14 @@ % % Algorithm description: % At time step k, define two matrix -% Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], +% 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 x(k-w+1), 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. % @@ -86,11 +88,12 @@ function update(obj, xold, yold, xnew, ynew) % Forget the oldest pair of snapshots (xold, yold), and remebers 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 - % Xw = [x(k-w+2),x(k-w+3),...,x(k+1)], Yw = [y(k-w+2),y(k-w+3),...,y(k+1)], + % 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 % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) + % Forget the oldest snapshot pair % compute Pk*xold matrix vector product beforehand Pkxold = obj.P*xold; % compute beta @@ -100,6 +103,7 @@ function update(obj, xold, yold, xnew, ynew) % compute Pko Pko = obj.P + (beta*Pkxold)*Pkxold'; + % Remember the newest snapshot pair % compute Pko*xnew matrix vector product beforehand Pkoxnew = Pko*xnew; % compute gamma diff --git a/matlab/WindowDMD.m~ b/matlab/WindowDMD.m~ deleted file mode 100644 index 3161379..0000000 --- a/matlab/WindowDMD.m~ +++ /dev/null @@ -1,128 +0,0 @@ -% WindowDMD is a class that implements window dynamic mode decomposition -% The time complexity (for one iteration) is O(n^2), and space complexity is -% O(n^2), where n is the state dimension -% -% Algorithm description: -% At time step k, define two matrix -% Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [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). -% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively -% by efficient rank-2 updating window DMD algorithm. -% -% Usage: -% wdmd = WindowDMD(n,w) -% wdmd.initialize(Xq,Yq) -% wdmd.update(xold, yold, xnew, ynew) -% [evals, modes] = wdmd.computemodes() -% -% properties: -% n: state dimension -% w: finite time window size -% timestep: number of snapshot pairs processed -% 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(Xq, Yq), initialize window DMD algorithm -% update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, -% and remeber new snapshot pair -% computemodes(), compute and return DMD eigenvalues and DMD mdoes -% -% 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", -% in production, 2017. To be submitted for publication, available on arXiv. -% -% Created: -% April 2017. -% -% To look up the documentation, type help WindowDMD - -classdef WindowDMD < handle - properties - n = 0; % state dimension - w = 0; % weighting factor - timestep = 0; % number of snapshots processed - 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) - % Creat an object for window DMD - % Usage: wdmd = WindowDMD(n,w) - if nargin == 2 - obj.n = n; - obj.w = w; - obj.A = zeros(n,n); - obj.P = zeros(n,n); - end - end - - function initialize(obj, Xq, Yq) - % Initialize WnlineDMD with q snapshot pairs stored in (Xq, Yq) - % Usage: wdmd.initialize(Xq,Yq) - q = length(Xq(1,:)); - if(obj.timestep == 0 && obj.w == q && obj.w >= obj.n+1) - obj.A = Yq*pinv(Xq); - obj.P = inv(Xq*Xq'); - end - obj.timestep = obj.timestep + q; - end - - function update(obj, xold, yold, xnew, ynew) - - % Update the DMD computation by sliding the finite time window forward - % Forget the oldest pair of snapshots (xold, yold), and remebers 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 - % Xw = [x(k-w+2),x(k-w+3),...,x(k+1)], Yw = [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 - % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) - % Usage: wdmd.update(xold, yold, xnew, ynew) - - % compute Pk*xold matrix vector product beforehand - Pxold = obj.P*xold; - % compute beta - beta = 1/(1-xold'*Pkxold); - % compute Pko - Pko = obj.P + (beta*Pkxold)*Pkxold'; - - - % compute M*xnew matrix vector product beforehand - Mxnew = obj.M*xnew; - % Compute gamma - gamma = 1/(1+xnew'*(Mxnew)); - % Compute Pk+1 - Pk1 = obj.M - gamma*(Mxnew*Mxnew'); - % compute P(k+1)*xold matrix vector product beforehand - Pk1xold = Pk1*xold; - % Compute beta - beta = 1/(1-xold'*(Pk1xold)); - - % Update A - obj.A = obj.B + gamma*((ynew-obj.B*xnew)*Mxnew'); - % Update B - obj.B = obj.A + beta*((-yold+obj.A*xold)*Pk1xold'); - % Update M - obj.M = Pk1 + beta*(Pk1xold*Pk1xold'); - - % time step + 1 - obj.timestep = obj.timestep + 1; - end - - function [evals, modes] = computemodes(obj) - % Compute and return DMD eigenvalues and DMD modes at current time step - % Usage: [evals, modes] = wdmd.computemodes() - [modes, evals] = eig(obj.A, 'vector'); - end - end -end \ No newline at end of file diff --git a/matlab/online_demo.m b/matlab/online_demo.m index b8f618d..b8c5437 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -5,16 +5,17 @@ % w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) % are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit % -% At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], +% 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, we would like to compute -% Ak = Yk*pinv(Xk). This can be done by brute-force batch DMD, +% Ak = Yk*pinv(Xk). At time step K+1, we need to include new snapshot pair x(k+1), y(k+1) +% This can be done by brute-force batch DMD, % and by efficient rank-1 updating online DMD algrithm. % % Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly % % Online DMD computes the DMD matrix by using efficient rank-1 update idea % -% We compare the performance of online DMD (with lambda=1,0.9) with the brute-force batch DMD +% We compare the performance of online DMD (with alpha=1,0.9) with the brute-force batch DMD % approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution % % Authors: @@ -74,7 +75,7 @@ elapsed_time = toc; fprintf('Batch DMD, elapsed time: %f seconds\n', elapsed_time) -% Online DMD lambda = 1 +% Online DMD alpha = 1 q = 20; evalsonlineDMD1 = zeros(n,m); % creat object and initialize with first q snapshot pairs @@ -87,9 +88,9 @@ evalsonlineDMD1(:,k) = log(eig(odmd.A))/dt; end elapsed_time = toc; -fprintf('Online DMD, lambda = 1, elapsed time: %f seconds\n', elapsed_time) +fprintf('Online DMD, alpha = 1, elapsed time: %f seconds\n', elapsed_time) -% Online DMD, lambda = 0.9 +% Online DMD, alpha = 0.9 q = 20; evalsonlineDMD09 = zeros(n,m); % creat object and initialize with first q snapshot pairs @@ -102,11 +103,11 @@ evalsonlineDMD09(:,k) = log(eig(odmd.A))/dt; end elapsed_time = toc; -fprintf('Online DMD, lambda = 0.9, elapsed time: %f seconds\n', elapsed_time) +fprintf('Online DMD, alpha = 0.9, elapsed time: %f seconds\n', elapsed_time) % visualize imaginary part of the continous time eigenvalues -% from true, batch, online (lambda=1), and online (lambda=0.9) +% from true, batch, online (alpha=1), and online (alpha=0.9) updateindex = q+1:m; figure, hold on plot(t,imag(evals(1,:)),'k-','LineWidth',2) @@ -115,7 +116,7 @@ plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') title('Imaginary part of eigenvalues','Interpreter','latex') -fl = legend('True','batch','online, $\lambda=1$','online, $\lambda=0.9$'); +fl = legend('True','batch','online, $\alpha=1$','online, $\alpha=0.9$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on diff --git a/matlab/window_demo.m b/matlab/window_demo.m index a824a2b..c9d611b 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -7,7 +7,10 @@ % % At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], % that contain the recent w snapshot pairs from a finite time window, -% we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, +% we would like to compute Ak = Yk*pinv(Xk). +% At time step k+1, we need to forget old snapshot pair x(k-w+1), y(k-w+1), +% and remember new snapshot pair xnew = x(k+1), ynew = y(k+1) +% This can be done by brute-force mini-batch DMD, % and by efficient rank-2 updating window DMD algrithm. % % Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index ae9ba8e..6e575f4 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -16,7 +16,7 @@ class OnlineDMD: by efficient rank-1 updating online DMD algrithm. Usage: - odmd = OnlineDMD(n,forgetting) + odmd = OnlineDMD(n,weighting) odmd.initialize(Xq,Yq) odmd.initilizeghost() odmd.update(x,y) @@ -24,7 +24,7 @@ class OnlineDMD: properties: n: state dimension - forgetting: forgetting factor between (0,1] + weighting: weighting factor between (0,1] timestep: number of snapshot pairs processed (i.e., the current time step) A: DMD matrix, size n by n P: Matrix that contains information about past snapshots, size n by n @@ -49,13 +49,13 @@ class OnlineDMD: To import the OnlineDMD class, add import online at head of Python scripts. To look up this documentation, type help(online.OnlineDMD) or online.OnlineDMD? """ - def __init__(self, n=0, forgetting=1, timestep=0, A=None, P=None): + def __init__(self, n=0, weighting=1, timestep=0, A=None, P=None): """ Creat an object for online DMD - Usage: odmd = OnlineDMD(n,forgetting) + Usage: odmd = OnlineDMD(n,weighting) """ self.n = n - self.forgetting = forgetting + self.weighting = weighting self.timestep = timestep if A is None or P is None: self.A = np.zeros([n,n]) @@ -71,13 +71,14 @@ def initialize(self, Xq, Yq): q = len(Xq[0,:]) Xqhat, Yqhat = np.zeros(Xq.shape), np.zeros(Yq.shape) if self.timestep == 0 and self.n <= q: - sqrtlambda = np.sqrt(self.forgetting) - # multiply forgetting factor with snapshots + #weight = np.sqrt(self.weighting)**range(q-1,-1,-1) + sqrtlambda = np.sqrt(self.weighting) + # multiply weighting factor with snapshots for i in range(q): Xqhat[:,i] = Xq[:,i]*sqrtlambda**(q-1-i) Yqhat[:,i] = Yq[:,i]*sqrtlambda**(q-1-i) self.A = Yqhat.dot(np.linalg.pinv(Xqhat)) - self.P = np.linalg.inv(Xqhat.dot(Xqhat.T))/self.forgetting + self.P = np.linalg.inv(Xqhat.dot(Xqhat.T))/self.weighting self.timestep += q def initializeghost(self): @@ -102,7 +103,7 @@ def update(self, x, y): # update A self.A += np.outer(y-self.A.dot(x),gamma*Px) # update P - self.P = (self.P - np.outer(gamma*Px,Px))/self.forgetting + self.P = (self.P - np.outer(gamma*Px,Px))/self.weighting # time step + 1 self.timestep += 1 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 8d3350c..6f4a4d9 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -14,6 +14,8 @@ class WindowDMD: 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 correponding to consecutive states z(k-1) and z(k). + At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), + and remember the 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 algrithm. @@ -84,6 +86,7 @@ def update(self, xold, yold, xnew, ynew): xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) Usage: wdmd.update(xold, yold, xnew, ynew) """ + # Forget the oldest snapshot pair # compute Pkxold matrix vector product beforehand Pkxold = self.P.dot(xold) # compute beta @@ -93,6 +96,7 @@ def update(self, xold, yold, xnew, ynew): # compute Pko Pko = self.P + np.outer(beta*Pkxold, Pkxold) + # Remember the newest snapshot pair # compute Pko*xnew matrix vector product beforehand Pkoxnew = Pko.dot(xnew) # compute gamma diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index b4041da..a3a2359 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -15,7 +15,7 @@ Online DMD computes the DMD matrix by using efficient rank-1 update idea -We compare the performance of online DMD (with lambda=1,0.9) with the brute-force batch DMD +We compare the performance of online DMD (with alpha=1,0.9) with the brute-force batch DMD approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution Authors: @@ -89,7 +89,7 @@ def dyn(x,t): print "Batch DMD, time = " + str(end-start) + " secs" -# Online DMD, forgetting = 1 +# Online DMD, weighting = 1 q = 20 evalsonlineDMD1 = np.empty((n,m),dtype=complex) odmd = dmdtools.OnlineDMD(n,1.0) @@ -99,10 +99,10 @@ def dyn(x,t): odmd.update(x[:,k],y[:,k]) evalsonlineDMD1[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt end = time.clock() -print "Online DMD, forgetting = 1, time = " + str(end-start) + " secs" +print "Online DMD, weighting = 1, time = " + str(end-start) + " secs" -# Online DMD, forgetting = 0.9 +# Online DMD, weighting = 0.9 q = 20 evalsonlineDMD09 = np.empty((n,m),dtype=complex) odmd = dmdtools.OnlineDMD(n,0.9) @@ -112,17 +112,17 @@ def dyn(x,t): odmd.update(x[:,k],y[:,k]) evalsonlineDMD09[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt end = time.clock() -print "Online DMD, forgetting = 0.9, time = " + str(end-start) + " secs" +print "Online DMD, weighting = 0.9, time = " + str(end-start) + " secs" -# visualize true, batch, online (forgettting=1,0.9) +# visualize true, batch, online (weighting=1,0.9) plt.figure() plt.rc('text', usetex=True) plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='batch',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, $\lambda$=1',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, $\lambda$=0.9',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, weighting=1',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, weighting=0.9',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im', fontsize=20) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index f78012a..8e58fcf 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -9,7 +9,9 @@ At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], that contain the recent w snapshot pairs from a finite time window, we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, -and by efficient rank-2 updating window DMD algrithm. +and by efficient rank-2 updating window DMD algrithm. +For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), +and remember the new snapshot pair xnew = x(k+1), ynew = y(k+1) Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly From 2ea2b4a142b983ed18e625e85218034e22730180 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 11:30:50 -0400 Subject: [PATCH 15/32] update python implementation and demo --- python/dmdtools/online.py | 8 ++------ python/scripts/online_demo.py | 6 +++--- python/scripts/window_demo.py | 4 ++-- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 6e575f4..c1cba78 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -71,12 +71,8 @@ def initialize(self, Xq, Yq): q = len(Xq[0,:]) Xqhat, Yqhat = np.zeros(Xq.shape), np.zeros(Yq.shape) if self.timestep == 0 and self.n <= q: - #weight = np.sqrt(self.weighting)**range(q-1,-1,-1) - sqrtlambda = np.sqrt(self.weighting) - # multiply weighting factor with snapshots - for i in range(q): - Xqhat[:,i] = Xq[:,i]*sqrtlambda**(q-1-i) - Yqhat[:,i] = Yq[:,i]*sqrtlambda**(q-1-i) + weight = np.sqrt(self.weighting)**range(q-1,-1,-1) + Xqhat, Yqhat = weight*Xq, weight*Yq self.A = Yqhat.dot(np.linalg.pinv(Xqhat)) self.P = np.linalg.inv(Xqhat.dot(Xqhat.T))/self.weighting self.timestep += q diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index a3a2359..af3a5fb 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -83,7 +83,7 @@ def dyn(x,t): evalsbatchDMD = np.empty((n,m),dtype=complex) start = time.clock() for k in range(q,m): - AbatchDMD[:,:,k] = y[:,:k].dot(np.linalg.pinv(x[:,:k])) + AbatchDMD[:,:,k] = y[:,:k+1].dot(np.linalg.pinv(x[:,:k+1])) evalsbatchDMD[:,k] = np.log(np.linalg.eigvals(AbatchDMD[:,:,k]))/dt end = time.clock() print "Batch DMD, time = " + str(end-start) + " secs" @@ -121,8 +121,8 @@ def dyn(x,t): plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='batch',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, weighting=1',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, weighting=0.9',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, wf=1',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, wf=0.9',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im', fontsize=20) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 8e58fcf..d77b416 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -85,7 +85,7 @@ def dyn(x,t): evalsminibatchDMD = np.empty((n,m),dtype=complex) start = time.clock() for k in range(w,m): - AminibatchDMD[:,:,k] = y[:,k-w:k].dot(np.linalg.pinv(x[:,k-w:k])) + AminibatchDMD[:,:,k] = y[:,k-w+1:k+1].dot(np.linalg.pinv(x[:,k-w+1:k+1])) evalsminibatchDMD[:,k] = np.log(np.linalg.eigvals(AminibatchDMD[:,:,k]))/dt end = time.clock() print "Mini-batch DMD, time = " + str(end-start) + " secs" @@ -104,7 +104,7 @@ def dyn(x,t): print "Window DMD, time = " + str(end-start) + " secs" -# visualize true, batch, window (forgettting=1,0.9) +# visualize true, batch, window plt.figure() plt.rc('text', usetex=True) plt.rc('font', family='serif') From 947214856cf54deeecbf5a89a8f3b7eccafea297 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 17:49:02 -0400 Subject: [PATCH 16/32] direct rank-2 update of window DMD --- matlab/OnlineDMD.m | 23 ++++++++++++----------- matlab/WindowDMD.m | 30 ++++++++++++------------------ matlab/online_demo.m | 20 +++++++++----------- matlab/window_demo.m | 11 +++++------ python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 35 +++++++++++++++-------------------- python/scripts/online_demo.py | 16 ++++++++-------- python/scripts/window_demo.py | 12 ++++++------ 8 files changed, 68 insertions(+), 81 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index dd998bb..8e27ba8 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -13,7 +13,7 @@ % by efficient rank-1 updating online DMD algorithm. % % Usage: -% odmd = OnlineDMD(n,alpha) +% odmd = OnlineDMD(n,weighting) % odmd.initialize(Xq,Yq) % odmd.initilizeghost() % odmd.update(x,y) @@ -21,7 +21,7 @@ % % properties: % n: state dimension -% alpha: weighting factor between 0 and 1 +% weighting: weighting factor between 0 and 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 @@ -42,7 +42,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. To be submitted for publication, available on arXiv. +% in production, 2017. Available on arXiv. % % Created: % April 2017. @@ -52,19 +52,19 @@ classdef OnlineDMD < handle properties n = 0; % state dimension - alpha = 1; % weighting factor + weighting = 1; % weighting factor timestep = 0; % number of snapshots processed A; % DMD matrix P; % matrix that contains information about past snapshots end methods - function obj = OnlineDMD(n,alpha) + function obj = OnlineDMD(n,weighting) % Creat an object for online DMD - % Usage: odmd = OnlineDMD(n,alpha) + % Usage: odmd = OnlineDMD(n,weighting) if nargin == 2 obj.n = n; - obj.alpha = alpha; + obj.weighting = weighting; obj.A = zeros(n,n); obj.P = zeros(n,n); end @@ -75,11 +75,11 @@ function initialize(obj, Xq, Yq) % Usage: odmd.initialize(Xq,Yq) q = length(Xq(1,:)); if(obj.timestep == 0 && q>=obj.n) - weight = (sqrt(obj.alpha)).^(q-1:-1:0); + 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.alpha; + obj.P = inv(Xq*Xq')/obj.weighting; end obj.timestep = obj.timestep + q; end @@ -87,8 +87,9 @@ function initialize(obj, Xq, Yq) 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/eps)*eye(obj.n); + obj.P = (1/epsilon)*eye(obj.n); end function update(obj, x, y) @@ -104,7 +105,7 @@ function update(obj, x, y) % Update A obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; % Update P - obj.P = (obj.P - (gamma*Px)*Px')/obj.alpha; + obj.P = (obj.P - (gamma*Px)*Px')/obj.weighting; % time step + 1 obj.timestep = obj.timestep + 1; end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index e81713d..8aef3b1 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -43,7 +43,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. To be submitted for publication, available on arXiv. +% in production, 2017. Available on arXiv. % % Created: % April 2017. @@ -93,25 +93,19 @@ function update(obj, xold, yold, xnew, ynew) % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) - % Forget the oldest snapshot pair - % compute Pk*xold matrix vector product beforehand - Pkxold = obj.P*xold; - % compute beta - beta = 1/(1-xold'*Pkxold); - % compute Ako - Ako = obj.A + (beta*(-yold+obj.A*xold))*Pkxold'; - % compute Pko - Pko = obj.P + (beta*Pkxold)*Pkxold'; - - % Remember the newest snapshot pair - % compute Pko*xnew matrix vector product beforehand - Pkoxnew = Pko*xnew; - % compute gamma - gamma = 1/(1+xnew'*Pkoxnew); + % direct rank-2 update + % define matrices + U = [xold, xnew]; V = [yold, ynew]; C = [-1,0;0,1]; + % compute PkU matrix vector product beforehand + PkU = obj.P*U; + % compute AkU matrix vector product beforehand + AkU = obj.A*U; + % compute Gamma + Gamma = inv(C+U'*PkU); % update A - obj.A = Ako + (gamma*(ynew - Ako*xnew))*Pkoxnew'; + obj.A = obj.A + (V-AkU)*(Gamma*PkU'); % update P - obj.P = Pko - (gamma*Pkoxnew)*Pkoxnew'; + obj.P = obj.P - PkU*(Gamma*PkU'); % time step + 1 obj.timestep = obj.timestep + 1; diff --git a/matlab/online_demo.m b/matlab/online_demo.m index b8c5437..ff1dc87 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -14,10 +14,9 @@ % Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly % % Online DMD computes the DMD matrix by using efficient rank-1 update idea -% % We compare the performance of online DMD (with alpha=1,0.9) with the brute-force batch DMD % approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -% +% Online DMD (weighting=1) and batch DMD should agree with each other (up to machine round-offer errors) % Authors: % Hao Zhang % Clarence W. Rowley @@ -25,7 +24,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. To be submitted for publication, available on arXiv. +% in production, 2017. Available on arXiv. % % Date created: April 2017 @@ -75,7 +74,7 @@ elapsed_time = toc; fprintf('Batch DMD, elapsed time: %f seconds\n', elapsed_time) -% Online DMD alpha = 1 +% Online DMD weighting = 1 q = 20; evalsonlineDMD1 = zeros(n,m); % creat object and initialize with first q snapshot pairs @@ -88,9 +87,9 @@ evalsonlineDMD1(:,k) = log(eig(odmd.A))/dt; end elapsed_time = toc; -fprintf('Online DMD, alpha = 1, elapsed time: %f seconds\n', elapsed_time) +fprintf('Online DMD, weighting = 1, elapsed time: %f seconds\n', elapsed_time) -% Online DMD, alpha = 0.9 +% Online DMD, weighting = 0.9 q = 20; evalsonlineDMD09 = zeros(n,m); % creat object and initialize with first q snapshot pairs @@ -103,20 +102,19 @@ evalsonlineDMD09(:,k) = log(eig(odmd.A))/dt; end elapsed_time = toc; -fprintf('Online DMD, alpha = 0.9, elapsed time: %f seconds\n', elapsed_time) +fprintf('Online DMD, weighting = 0.9, elapsed time: %f seconds\n', elapsed_time) % visualize imaginary part of the continous time eigenvalues -% from true, batch, online (alpha=1), and online (alpha=0.9) +% from true, batch, online (rho=1), and online (rho=0.9) updateindex = q+1:m; figure, hold on plot(t,imag(evals(1,:)),'k-','LineWidth',2) plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) -xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') -title('Imaginary part of eigenvalues','Interpreter','latex') -fl = legend('True','batch','online, $\alpha=1$','online, $\alpha=0.9$'); +xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') +fl = legend('True','Batch','Online, $wf=1$','Online, $wf=0.9$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on diff --git a/matlab/window_demo.m b/matlab/window_demo.m index c9d611b..b6324e2 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -19,7 +19,8 @@ % % We compare the performance of window DMD with the brute-force mini-batch DMD % approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -% +% They should agree with each other (up to machine round-offer errors) +% % Authors: % Hao Zhang % Clarence W. Rowley @@ -27,7 +28,7 @@ % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, % ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. To be submitted for publication, available on arXiv. +% in production, 2017. Available on arXiv. % % Date created: April 2017 @@ -80,7 +81,6 @@ % window DMD -w = 20; evalswindowDMD = zeros(n,m); % creat object and initialize with first w snapshot pairs wdmd = WindowDMD(n,w); @@ -102,9 +102,8 @@ plot(t,imag(evals(1,:)),'k-','LineWidth',2) plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) -xlabel('Time','Interpreter','latex'), ylabel('Im','Interpreter','latex') -title('Imaginary part of eigenvalues','Interpreter','latex') -fl = legend('True','mini-batch','window'); +xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') +fl = legend('True','Mini-batch, $w=20$','Window, $w=20$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index c1cba78..2e55e06 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -42,7 +42,7 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. To be submitted for publication, available on arXiv. + in production, 2017. Available on arXiv. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 6f4a4d9..4e211d6 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -44,7 +44,7 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. To be submitted for publication, available on arXiv. + in production, 2017. Available on arXiv. Date created: April 2017 @@ -86,25 +86,20 @@ def update(self, xold, yold, xnew, ynew): xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) Usage: wdmd.update(xold, yold, xnew, ynew) """ - # Forget the oldest snapshot pair - # compute Pkxold matrix vector product beforehand - Pkxold = self.P.dot(xold) - # compute beta - beta = 1.0/(1-xold.T.dot(Pkxold)) - # compute Ako - Ako = self.A + np.outer(beta*(-yold+self.A.dot(xold)),Pkxold) - # compute Pko - Pko = self.P + np.outer(beta*Pkxold, Pkxold) - - # Remember the newest snapshot pair - # compute Pko*xnew matrix vector product beforehand - Pkoxnew = Pko.dot(xnew) - # compute gamma - gamma = 1.0/(1+xnew.T.dot(Pkoxnew)) - # update Ak - self.A = Ako + np.outer(gamma*(ynew-Ako.dot(xnew)),Pkoxnew) - # update Pk - self.P = Pko - np.outer(gamma*Pkoxnew, Pkoxnew) + # direct rank-2 update + # define matrices + U, V = np.column_stack((xold, xnew)), np.column_stack((yold, ynew)) + C = np.diag([-1,1]) + # compute PkU matrix vector product beforehand + PkU = self.P.dot(U) + # compute AkU matrix vector product beforehand + AkU = self.A.dot(U) + # compute Gamma + Gamma = np.linalg.inv(C+U.T.dot(PkU)) + # update A + self.A += (V-AkU).dot(Gamma).dot(PkU.T) + # update P + self.P += -PkU.dot(Gamma).dot(PkU.T) # time step + 1 self.timestep += 1 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index af3a5fb..1707b9a 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -15,8 +15,9 @@ Online DMD computes the DMD matrix by using efficient rank-1 update idea -We compare the performance of online DMD (with alpha=1,0.9) with the brute-force batch DMD +We compare the performance of online DMD (with weighting=1,0.9) with the brute-force batch DMD approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution +Online DMD (weighting=1) and batch DMD should agree with each other (up to machine round-offer errors) Authors: Hao Zhang @@ -25,7 +26,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. To be submitted for publication, available on arXiv. + in production, 2017. Available on arXiv. Date created: April 2017 """ @@ -119,14 +120,13 @@ def dyn(x,t): plt.figure() plt.rc('text', usetex=True) plt.rc('font', family='serif') -plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='batch',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='online, wf=1',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='online, wf=0.9',linewidth=2.0) +plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='Batch',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='Online, wf=1',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='Online, wf=0.9',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) -plt.ylabel('Im', fontsize=20) -plt.title('Imignary part of eigenvalues', fontsize=20) +plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) plt.legend(loc='best', fontsize=20, shadow=True) plt.xlim([0,10]) plt.ylim([1,2]) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index d77b416..528f83b 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -19,6 +19,7 @@ We compare the performance of window DMD with the brute-force mini-batch DMD approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution +They should agree with each other (up to machine round-offer errors) Authors: Hao Zhang @@ -27,7 +28,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. To be submitted for publication, available on arXiv. + in production, 2017. Available on arXiv. Date created: April 2017 """ @@ -108,13 +109,12 @@ def dyn(x,t): plt.figure() plt.rc('text', usetex=True) plt.rc('font', family='serif') -plt.plot(t, np.imag(evals[0,:]), 'k-',label='true',linewidth=2.0) -plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='mini-batch, w=20',linewidth=2.0) -plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='window, w=20',linewidth=2.0) +plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) +plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, w=20',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, w=20',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) -plt.ylabel('Im', fontsize=20) -plt.title('Imignary part of eigenvalues', fontsize=20) +plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) plt.legend(loc='best', fontsize=20, shadow=True) plt.xlim([0,10]) plt.ylim([1,2]) From 36078538580b4ec604d584796d584d3ddab5c5d4 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 17:56:26 -0400 Subject: [PATCH 17/32] update plots --- python/scripts/online_demo.py | 4 ++-- python/scripts/window_demo.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 1707b9a..186a67b 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -122,8 +122,8 @@ def dyn(x,t): plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='Batch',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='Online, wf=1',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='Online, wf=0.9',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='Online, $wf=1$',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='Online, $wf=0.9$',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 528f83b..d17e268 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -110,8 +110,8 @@ def dyn(x,t): plt.rc('text', usetex=True) plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) -plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, w=20',linewidth=2.0) -plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, w=20',linewidth=2.0) +plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, $w=20$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, $w=20$',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) From a1a8727fc6051fa6446acefef0a6c3732f5bd252 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 20:02:16 -0400 Subject: [PATCH 18/32] fix documentation matrix vector definition --- matlab/WindowDMD.m | 6 +++--- matlab/window_demo.m | 4 ++-- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 6 +++--- python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 4 ++-- 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 8aef3b1..6d33f47 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -10,7 +10,7 @@ % 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 x(k-w+1), y(k-w+1), +% At time step k+1, we need to forget old snapshot pair x(k-w), y(k-w), % 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. @@ -33,7 +33,7 @@ % update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, % and remeber new snapshot pair % At time k+1, if 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 xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) +% we should take xold = x(k-w), yold = y(k-w), xnew = x(k+1), ynew = y(k+1) % computemodes(), compute and return DMD eigenvalues and DMD mdoes % % Authors: @@ -90,7 +90,7 @@ function update(obj, xold, yold, xnew, ynew) % 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 - % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) + % xold = x(k-w), yold = y(k-w), xnew = x(k), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) % direct rank-2 update diff --git a/matlab/window_demo.m b/matlab/window_demo.m index b6324e2..de4bad3 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -5,10 +5,10 @@ % w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) % are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit % -% At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], +% 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, % we would like to compute Ak = Yk*pinv(Xk). -% At time step k+1, we need to forget old snapshot pair x(k-w+1), y(k-w+1), +% At time step k+1, we need to forget old snapshot pair x(k-w), y(k-w), % and remember new snapshot pair xnew = x(k+1), ynew = y(k+1) % This can be done by brute-force mini-batch DMD, % and by efficient rank-2 updating window DMD algrithm. diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 2e55e06..dca9ad1 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -7,7 +7,7 @@ class OnlineDMD: O(n^2), where n is the state dimension. Algorithm description: - At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], + 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) diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 4e211d6..43a1be3 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -8,13 +8,13 @@ class WindowDMD: Algorithm description: At time step k, define two matrix - Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], + 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 correponding to consecutive states z(k-1) and z(k). - At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), + At time k+1, we need to forget the old snapshot pair xold = x(k-w), yold = y(k-w), and remember the 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 algrithm. @@ -83,7 +83,7 @@ def update(self, xold, yold, xnew, ynew): 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 - xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) + xold = x(k-w), yold = y(k-w), xnew = x(k+1), ynew = y(k+1) Usage: wdmd.update(xold, yold, xnew, ynew) """ # direct rank-2 update diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 186a67b..b01e29d 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -6,7 +6,7 @@ w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit -At time step k, define two matrix Xk = [x(1),x(2),...,x(k)], Yk = [y(1),y(2),...,y(k)], +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, we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force batch DMD, and by efficient rank-1 updating online DMD algrithm. diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index d17e268..e21076f 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -6,11 +6,11 @@ w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit -At time step k, define two matrix Xk = [x(k-w+1),x(k-w+2),...,x(k)], Yk = [y(k-w+1),y(k-w+2),...,y(k)], +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, we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, and by efficient rank-2 updating window DMD algrithm. -For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), +For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w), yold = y(k-w), and remember the new snapshot pair xnew = x(k+1), ynew = y(k+1) Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly From b8e351ba4b0597277c40240224cf70a190b29805 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 20:18:30 -0400 Subject: [PATCH 19/32] fix window DMD documentation --- matlab/WindowDMD.m | 6 +++--- matlab/window_demo.m | 2 +- python/dmdtools/window.py | 4 ++-- python/scripts/window_demo.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 6d33f47..9ae8694 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -10,7 +10,7 @@ % 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 x(k-w), y(k-w), +% 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. @@ -33,7 +33,7 @@ % update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, % and remeber new snapshot pair % At time k+1, if 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 xold = x(k-w), yold = y(k-w), xnew = x(k+1), ynew = y(k+1) +% we should take xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % computemodes(), compute and return DMD eigenvalues and DMD mdoes % % Authors: @@ -90,7 +90,7 @@ function update(obj, xold, yold, xnew, ynew) % 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 - % xold = x(k-w), yold = y(k-w), xnew = x(k), ynew = y(k+1) + % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) % Usage: wdmd.update(xold, yold, xnew, ynew) % direct rank-2 update diff --git a/matlab/window_demo.m b/matlab/window_demo.m index de4bad3..ffdcf13 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -8,7 +8,7 @@ % 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, % we would like to compute Ak = Yk*pinv(Xk). -% At time step k+1, we need to forget old snapshot pair x(k-w), y(k-w), +% 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) % This can be done by brute-force mini-batch DMD, % and by efficient rank-2 updating window DMD algrithm. diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 43a1be3..3b5b636 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -14,7 +14,7 @@ class WindowDMD: 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 correponding to consecutive states z(k-1) and z(k). - At time k+1, we need to forget the old snapshot pair xold = x(k-w), yold = y(k-w), + At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), and remember the 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 algrithm. @@ -83,7 +83,7 @@ def update(self, xold, yold, xnew, ynew): 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 - xold = x(k-w), yold = y(k-w), xnew = x(k+1), ynew = y(k+1) + xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) Usage: wdmd.update(xold, yold, xnew, ynew) """ # direct rank-2 update diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index e21076f..2d475e3 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -10,7 +10,7 @@ that contain the recent w snapshot pairs from a finite time window, we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, and by efficient rank-2 updating window DMD algrithm. -For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w), yold = y(k-w), +For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), and remember the new snapshot pair xnew = x(k+1), ynew = y(k+1) Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly From b3bfbb19e0c5237180d25fd1d3f5e67cb917ac86 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 21:35:15 -0400 Subject: [PATCH 20/32] fix demo ode solver --- matlab/online_demo.m | 16 ++++++---------- matlab/window_demo.m | 16 +++++++--------- python/scripts/online_demo.py | 5 ++--- python/scripts/window_demo.py | 12 ++++++------ 4 files changed, 21 insertions(+), 28 deletions(-) diff --git a/matlab/online_demo.m b/matlab/online_demo.m index ff1dc87..fdb0866 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -32,15 +32,13 @@ epsilon = 1e-1; dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x; % generate data -tspan = [0 10]; -x0 = [1;0]; -[t,x] = ode45(dyn, tspan, x0); -% interpolate uniform time step dt = 1e-1; -time = 0:dt:max(tspan); -xq = interp1(t,x,time); xq = xq'; +tspan = 0:dt:10; +x0 = [1;0]; +[tq,xq] = ode45(dyn, tspan, x0); % extract snapshot pairs -x = xq(:,1:end-1); y = xq(:,2:end); t = time(2:end); +xq = xq'; tq = tq'; +x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end); % true dynamics, eigenvalues [n, m] = size(x); A = zeros(n,n,m); @@ -63,7 +61,7 @@ % batch DMD -q = 20; +q = 10; AbatchDMD = zeros(n,n,m); evalsbatchDMD = zeros(n,m); tic @@ -75,7 +73,6 @@ fprintf('Batch DMD, elapsed time: %f seconds\n', elapsed_time) % Online DMD weighting = 1 -q = 20; evalsonlineDMD1 = zeros(n,m); % creat object and initialize with first q snapshot pairs odmd = OnlineDMD(n,1); @@ -90,7 +87,6 @@ fprintf('Online DMD, weighting = 1, elapsed time: %f seconds\n', elapsed_time) % Online DMD, weighting = 0.9 -q = 20; evalsonlineDMD09 = zeros(n,m); % creat object and initialize with first q snapshot pairs odmd = OnlineDMD(n,0.9); diff --git a/matlab/window_demo.m b/matlab/window_demo.m index ffdcf13..73aedf5 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -36,15 +36,13 @@ epsilon = 1e-1; dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x; % generate data -tspan = [0 10]; -x0 = [1;0]; -[t,x] = ode45(dyn, tspan, x0); -% interpolate uniform time step dt = 1e-1; -time = 0:dt:max(tspan); -xq = interp1(t,x,time); xq = xq'; +tspan = 0:dt:10; +x0 = [1;0]; +[tq,xq] = ode45(dyn, tspan, x0); % extract snapshot pairs -x = xq(:,1:end-1); y = xq(:,2:end); t = time(2:end); +xq = xq'; tq = tq'; +x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end); % true dynamics, eigenvalues [n, m] = size(x); A = zeros(n,n,m); @@ -67,7 +65,7 @@ % mini-batch DMD -w = 20; % storage time window size, store recent w snapshot pairs +w = 10; % storage time window size, store recent w snapshot pairs AminibatchDMD = zeros(n,n,m); evalsminibatchDMD = zeros(n,m); % mini-batch DMD @@ -103,7 +101,7 @@ plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') -fl = legend('True','Mini-batch, $w=20$','Window, $w=20$'); +fl = legend('True','Mini-batch, $w=10$','Window, $w=10$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index b01e29d..cf7b114 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -53,6 +53,7 @@ def dyn(x,t): dt = 0.1 x0 = [1,0] xsol = odeint(dyn,x0,tspan).T +# extract snapshots x, y = xsol[:,:-1], xsol[:,1:] t = tspan[:-1] # true dynamics, true eigenvalues @@ -79,7 +80,7 @@ def dyn(x,t): # batch DMD -q = 20 +q = 10 AbatchDMD = np.empty((n,n,m)) evalsbatchDMD = np.empty((n,m),dtype=complex) start = time.clock() @@ -91,7 +92,6 @@ def dyn(x,t): # Online DMD, weighting = 1 -q = 20 evalsonlineDMD1 = np.empty((n,m),dtype=complex) odmd = dmdtools.OnlineDMD(n,1.0) odmd.initialize(x[:,:q],y[:,:q]) @@ -104,7 +104,6 @@ def dyn(x,t): # Online DMD, weighting = 0.9 -q = 20 evalsonlineDMD09 = np.empty((n,m),dtype=complex) odmd = dmdtools.OnlineDMD(n,0.9) odmd.initialize(x[:,:q],y[:,:q]) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 2d475e3..525490d 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -55,6 +55,7 @@ def dyn(x,t): dt = 0.1 x0 = [1,0] xsol = odeint(dyn,x0,tspan).T +# extract snapshots x, y = xsol[:,:-1], xsol[:,1:] t = tspan[:-1] # true dynamics, true eigenvalues @@ -80,8 +81,8 @@ def dyn(x,t): plt.show() -# mini-batch DMD, w = 20 -w = 20 +# mini-batch DMD, w = 10 +w = 10 AminibatchDMD = np.empty((n,n,m)) evalsminibatchDMD = np.empty((n,m),dtype=complex) start = time.clock() @@ -92,8 +93,7 @@ def dyn(x,t): print "Mini-batch DMD, time = " + str(end-start) + " secs" -# Window DMD, w = 20 -w = 20 +# Window DMD, w = 10 evalswindowDMD = np.empty((n,m),dtype=complex) wdmd = dmdtools.WindowDMD(n,w) wdmd.initialize(x[:,:w],y[:,:w]) @@ -110,8 +110,8 @@ def dyn(x,t): plt.rc('text', usetex=True) plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) -plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, $w=20$',linewidth=2.0) -plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, $w=20$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, $w=10$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, $w=10$',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) From 14aa0d7abcd2020097b18202f999ad465e5bab1f Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 21:53:54 -0400 Subject: [PATCH 21/32] update demo --- matlab/online_demo.m | 14 +++++++------- matlab/window_demo.m | 12 ++++++------ python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 2 +- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/matlab/online_demo.m b/matlab/online_demo.m index fdb0866..b7524af 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -38,20 +38,20 @@ [tq,xq] = ode45(dyn, tspan, x0); % extract snapshot pairs xq = xq'; tq = tq'; -x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end); +x = xq(:,1:end-1); y = xq(:,2:end); time = tq(2:end); % true dynamics, eigenvalues [n, m] = size(x); A = zeros(n,n,m); evals = zeros(n,m); for k = 1:m - A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics + A(:,:,k) = [0, 1+epsilon*time(k); -(1+epsilon*time(k)),0]; % continuous time dynamics evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues end % visualize snapshots figure, hold on -plot(time,xq(1,:),'x-',time,xq(2,:),'o-','LineWidth',2) +plot(tq,xq(1,:),'x-',tq,xq(2,:),'o-','LineWidth',2) xlabel('Time','Interpreter','latex') title('Snapshots','Interpreter','latex') fl = legend('$x_1(t)$','$x_2(t)$'); @@ -105,10 +105,10 @@ % from true, batch, online (rho=1), and online (rho=0.9) updateindex = q+1:m; figure, hold on -plot(t,imag(evals(1,:)),'k-','LineWidth',2) -plot(t(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) -plot(t(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) -plot(t(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) +plot(time,imag(evals(1,:)),'k-','LineWidth',2) +plot(time(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) +plot(time(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) +plot(time(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') fl = legend('True','Batch','Online, $wf=1$','Online, $wf=0.9$'); set(fl,'Interpreter','latex','Location','northwest'); diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 73aedf5..8dec9a7 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -42,20 +42,20 @@ [tq,xq] = ode45(dyn, tspan, x0); % extract snapshot pairs xq = xq'; tq = tq'; -x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end); +x = xq(:,1:end-1); y = xq(:,2:end); time = tq(2:end); % true dynamics, eigenvalues [n, m] = size(x); A = zeros(n,n,m); evals = zeros(n,m); for k = 1:m - A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics + A(:,:,k) = [0, 1+epsilon*time(k); -(1+epsilon*time(k)),0]; % continuous time dynamics evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues end % visualize snapshots figure, hold on -plot(time,xq(1,:),'x-',time,xq(2,:),'o-','LineWidth',2) +plot(tq,xq(1,:),'x-',tq,xq(2,:),'o-','LineWidth',2) xlabel('Time','Interpreter','latex') title('Snapshots','Interpreter','latex') fl = legend('$x_1(t)$','$x_2(t)$'); @@ -97,9 +97,9 @@ % from true, mini-batch, and window updateindex = w+1:m; figure, hold on -plot(t,imag(evals(1,:)),'k-','LineWidth',2) -plot(t(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) -plot(t(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) +plot(time,imag(evals(1,:)),'k-','LineWidth',2) +plot(time(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) +plot(time(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') fl = legend('True','Mini-batch, $w=10$','Window, $w=10$'); set(fl,'Interpreter','latex','Location','northwest'); diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index cf7b114..7812dd0 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -55,7 +55,7 @@ def dyn(x,t): xsol = odeint(dyn,x0,tspan).T # extract snapshots x, y = xsol[:,:-1], xsol[:,1:] -t = tspan[:-1] +t = tspan[1:] # true dynamics, true eigenvalues n, m = len(x[:,0]), len(x[0,:]) A = np.empty((n,n,m)) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 525490d..8b4dee4 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -57,7 +57,7 @@ def dyn(x,t): xsol = odeint(dyn,x0,tspan).T # extract snapshots x, y = xsol[:,:-1], xsol[:,1:] -t = tspan[:-1] +t = tspan[1:] # true dynamics, true eigenvalues n, m = len(x[:,0]), len(x[0,:]) A = np.empty((n,n,m)) From ccd71f11d77520e2aef6f25bcbed1b8f36328419 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 22:36:32 -0400 Subject: [PATCH 22/32] fix online small weighting --- matlab/OnlineDMD.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 8e27ba8..f126de1 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -105,7 +105,7 @@ function update(obj, x, y) % Update A obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; % Update P - obj.P = (obj.P - (gamma*Px)*Px')/obj.weighting; + obj.P = (obj.P - gamma*(Px*Px'))/obj.weighting; % time step + 1 obj.timestep = obj.timestep + 1; end From 2d1691d7f9482be1e39073ac88b4021781496a46 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Fri, 19 May 2017 22:49:49 -0400 Subject: [PATCH 23/32] fix python online positive definite matrix computation --- matlab/OnlineDMD.m | 2 +- python/dmdtools/online.py | 6 +++--- python/scripts/window_demo.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index f126de1..ebf49fc 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -104,7 +104,7 @@ function update(obj, x, y) gamma = 1/(1+x'*Px); % Update A obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; - % Update P + % Update P, Px*Px' to ensure positive definite obj.P = (obj.P - gamma*(Px*Px'))/obj.weighting; % time step + 1 obj.timestep = obj.timestep + 1; diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index dca9ad1..b39ad53 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -97,9 +97,9 @@ def update(self, x, y): # compute gamma gamma = 1.0/(1 + x.T.dot(Px)) # update A - self.A += np.outer(y-self.A.dot(x),gamma*Px) - # update P - self.P = (self.P - np.outer(gamma*Px,Px))/self.weighting + self.A += np.outer(gamma*(y-self.A.dot(x)),Px) + # update P, Px*Px' to ensure positive definite + self.P = (self.P - gamma*np.outer(Px,Px))/self.weighting # time step + 1 self.timestep += 1 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 8b4dee4..69fc7b3 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -82,7 +82,7 @@ def dyn(x,t): # mini-batch DMD, w = 10 -w = 10 +w = 3 AminibatchDMD = np.empty((n,n,m)) evalsminibatchDMD = np.empty((n,m),dtype=complex) start = time.clock() From d00d7f23e4ae935ef05208f0f484ae5651a4f6fa Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Sun, 21 May 2017 13:58:31 -0400 Subject: [PATCH 24/32] fix online, ensure symmetric positive definite Pk --- matlab/OnlineDMD.m | 4 +++- matlab/WindowDMD.m | 2 +- python/dmdtools/online.py | 4 +++- python/scripts/window_demo.py | 2 +- 4 files changed, 8 insertions(+), 4 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index ebf49fc..e849e6b 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -104,8 +104,10 @@ function update(obj, x, y) gamma = 1/(1+x'*Px); % Update A obj.A = obj.A + (gamma*(y-obj.A*x))*Px'; - % Update P, Px*Px' to ensure positive definite + % 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 diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 9ae8694..0d0ed6c 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -95,7 +95,7 @@ function update(obj, xold, yold, xnew, ynew) % direct rank-2 update % define matrices - U = [xold, xnew]; V = [yold, ynew]; C = [-1,0;0,1]; + U = [xold, xnew]; V = [yold, ynew]; C = diag([-1,1]); % compute PkU matrix vector product beforehand PkU = obj.P*U; % compute AkU matrix vector product beforehand diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index b39ad53..b1227ba 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -98,8 +98,10 @@ def update(self, x, y): gamma = 1.0/(1 + x.T.dot(Px)) # update A self.A += np.outer(gamma*(y-self.A.dot(x)),Px) - # update P, Px*Px' to ensure positive definite + # update P, group Px*Px' to ensure positive definite self.P = (self.P - gamma*np.outer(Px,Px))/self.weighting + # ensure P is SPD by taking its symmetric part + self.P = (self.P + self.P.T)/2 # time step + 1 self.timestep += 1 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 69fc7b3..8b4dee4 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -82,7 +82,7 @@ def dyn(x,t): # mini-batch DMD, w = 10 -w = 3 +w = 10 AminibatchDMD = np.empty((n,n,m)) evalsminibatchDMD = np.empty((n,m),dtype=complex) start = time.clock() From 8f1f403bb04f64e6cfcf0d09e95fbbc677af43e9 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Sun, 21 May 2017 14:03:59 -0400 Subject: [PATCH 25/32] fix window, ensure symmetric positive definite Pk --- matlab/WindowDMD.m | 2 ++ python/dmdtools/window.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 0d0ed6c..fc8751f 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -106,6 +106,8 @@ function update(obj, xold, yold, xnew, ynew) obj.A = obj.A + (V-AkU)*(Gamma*PkU'); % update P obj.P = obj.P - PkU*(Gamma*PkU'); + % ensure P is SPD by taking its symmetric part + obj.P = (obj.P+(obj.P)')/2; % time step + 1 obj.timestep = obj.timestep + 1; diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 3b5b636..d889614 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -100,6 +100,8 @@ def update(self, xold, yold, xnew, ynew): self.A += (V-AkU).dot(Gamma).dot(PkU.T) # update P self.P += -PkU.dot(Gamma).dot(PkU.T) + # ensure P is SPD by taking its symmetric part + self.P = (self.P + self.P.T)/2 # time step + 1 self.timestep += 1 From d653c6f7875bd4debdcc959dfa3da34e38a94862 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Tue, 23 May 2017 14:43:19 -0400 Subject: [PATCH 26/32] add weighted version for window DMD --- README.md | 8 +++---- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 41 ++++++++++++++++++++--------------- matlab/online_demo.m | 6 ++--- matlab/window_demo.m | 30 ++++++++++++++++++------- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 35 +++++++++++++++++------------- python/scripts/online_demo.py | 6 ++--- python/scripts/window_demo.py | 27 +++++++++++++++++------ 9 files changed, 98 insertions(+), 59 deletions(-) diff --git a/README.md b/README.md index 1031a28..f0d1fa7 100644 --- a/README.md +++ b/README.md @@ -4,12 +4,12 @@ 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 online DMD) -4. Window DMD +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 online DMD) -5. Window DMD \ No newline at end of file +4. Online DMD (including weighted version) +5. Window DMD (including weighted version) \ No newline at end of file diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index e849e6b..59286e5 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -74,7 +74,7 @@ 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 && q>=obj.n) + if(obj.timestep == 0 && rank(Xq) == obj.n) weight = (sqrt(obj.weighting)).^(q-1:-1:0); Xq = Xq.*weight; Yq = Yq.*weight; diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index fc8751f..939b82b 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -16,20 +16,21 @@ % by efficient rank-2 updating window DMD algorithm. % % Usage: -% wdmd = WindowDMD(n,w) -% wdmd.initialize(Xq,Yq) +% wdmd = WindowDMD(n,w,weighting) +% wdmd.initialize(Xw,Yw) % wdmd.update(xold, yold, xnew, ynew) % [evals, modes] = wdmd.computemodes() % % properties: % n: state dimension % w: finite time window size +% weighting: weighting factor % timestep: number of snapshot pairs processed % 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(Xq, Yq), initialize window DMD algorithm +% initialize(Xw, Yw), initialize window DMD algorithm % update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, % and remeber new snapshot pair % At time k+1, if 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)], @@ -53,31 +54,36 @@ classdef WindowDMD < handle properties n = 0; % state dimension - w = 0; % weighting factor + w = 0; % window size + weighting = 1; % weighting factor timestep = 0; % number of snapshots processed 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) + function obj = WindowDMD(n, w, weighting) % Creat an object for window DMD - % Usage: wdmd = WindowDMD(n,w) - if nargin == 2 + % Usage: wdmd = WindowDMD(n,w,weighting) + if nargin == 3 obj.n = n; obj.w = w; + obj.weighting = weighting; obj.A = zeros(n,n); obj.P = zeros(n,n); end end - function initialize(obj, Xq, Yq) - % Initialize WnlineDMD with q snapshot pairs stored in (Xq, Yq) - % Usage: wdmd.initialize(Xq,Yq) - q = length(Xq(1,:)); - if(obj.timestep == 0 && obj.w == q && obj.w >= obj.n+1) - obj.A = Yq*pinv(Xq); - obj.P = inv(Xq*Xq'); + function initialize(obj, Xw, Yw) + % Initialize WnlineDMD with w snapshot pairs stored in (Xw, Yw) + % Usage: wdmd.initialize(Xw,Yw) + 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 @@ -95,17 +101,18 @@ function update(obj, xold, yold, xnew, ynew) % direct rank-2 update % define matrices - U = [xold, xnew]; V = [yold, ynew]; C = diag([-1,1]); + U = [xold, xnew]; V = [yold, ynew]; + C = diag([-(obj.weighting)^(obj.w),1]); % compute PkU matrix vector product beforehand PkU = obj.P*U; % compute AkU matrix vector product beforehand AkU = obj.A*U; % compute Gamma - Gamma = inv(C+U'*PkU); + 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.P = (obj.P - PkU*(Gamma*PkU'))/obj.weighting; % ensure P is SPD by taking its symmetric part obj.P = (obj.P+(obj.P)')/2; diff --git a/matlab/online_demo.m b/matlab/online_demo.m index b7524af..b2f03c0 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -87,7 +87,7 @@ fprintf('Online DMD, weighting = 1, elapsed time: %f seconds\n', elapsed_time) % Online DMD, weighting = 0.9 -evalsonlineDMD09 = zeros(n,m); +evalsonlineDMD2 = zeros(n,m); % creat object and initialize with first q snapshot pairs odmd = OnlineDMD(n,0.9); odmd.initialize(x(:,1:q),y(:,1:q)); @@ -95,7 +95,7 @@ tic for k = q+1:m odmd.update(x(:,k),y(:,k)); - evalsonlineDMD09(:,k) = log(eig(odmd.A))/dt; + evalsonlineDMD2(:,k) = log(eig(odmd.A))/dt; end elapsed_time = toc; fprintf('Online DMD, weighting = 0.9, elapsed time: %f seconds\n', elapsed_time) @@ -108,7 +108,7 @@ plot(time,imag(evals(1,:)),'k-','LineWidth',2) plot(time(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) plot(time(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) -plot(time(updateindex),imag(evalsonlineDMD09(1,updateindex)),'-','LineWidth',2) +plot(time(updateindex),imag(evalsonlineDMD2(1,updateindex)),'--','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') fl = legend('True','Batch','Online, $wf=1$','Online, $wf=0.9$'); set(fl,'Interpreter','latex','Location','northwest'); diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 8dec9a7..f9967f6 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -75,23 +75,36 @@ evalsminibatchDMD(:,k) = log(eig(AminibatchDMD(:,:,k)))/dt; end elapsed_time = toc; -fprintf('Mini-batch DMD, elapsed time: %f seconds\n', elapsed_time) +fprintf('Mini-batch DMD, w=10, elapsed time: %f seconds\n', elapsed_time) -% window DMD -evalswindowDMD = zeros(n,m); +% window DMD, weighting = 1 +evalswindowDMD1 = zeros(n,m); % creat object and initialize with first w snapshot pairs -wdmd = WindowDMD(n,w); +wdmd = WindowDMD(n,w,1); wdmd.initialize(x(:,1:w), y(:,1:w)); % window DMD tic for k = w+1:m wdmd.update(x(:,k-w), y(:,k-w), x(:,k), y(:,k)); - evalswindowDMD(:,k) = log(eig(wdmd.A))/dt; + evalswindowDMD1(:,k) = log(eig(wdmd.A))/dt; end elapsed_time = toc; -fprintf('Window DMD, elapsed time: %f seconds\n', elapsed_time) +fprintf('Window DMD, w = 10, weighting = 1, elapsed time: %f seconds\n', elapsed_time) +% window DMD, weighting = 0.5 +evalswindowDMD2 = zeros(n,m); +% creat object and initialize with first w snapshot pairs +wdmd = WindowDMD(n,w,0.5); +wdmd.initialize(x(:,1:w), y(:,1:w)); +% window DMD +tic +for k = w+1:m + wdmd.update(x(:,k-w), y(:,k-w), x(:,k), y(:,k)); + evalswindowDMD2(:,k) = log(eig(wdmd.A))/dt; +end +elapsed_time = toc; +fprintf('Window DMD, w = 10, weighting = 0.5, elapsed time: %f seconds\n', elapsed_time) % visualize imaginary part of the continous time eigenvalues % from true, mini-batch, and window @@ -99,9 +112,10 @@ figure, hold on plot(time,imag(evals(1,:)),'k-','LineWidth',2) plot(time(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) -plot(time(updateindex),imag(evalswindowDMD(1,updateindex)),'--','LineWidth',2) +plot(time(updateindex),imag(evalswindowDMD1(1,updateindex)),'--','LineWidth',2) +plot(time(updateindex),imag(evalswindowDMD2(1,updateindex)),'--','LineWidth',2) xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') -fl = legend('True','Mini-batch, $w=10$','Window, $w=10$'); +fl = legend('True','Mini-batch, $w=10$','Window, $w=10$, $wf=1$','Window, $w=10$, $wf=0.5$'); set(fl,'Interpreter','latex','Location','northwest'); ylim([1,2]), xlim([0,10]) box on diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index b1227ba..4dfb5bc 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -70,7 +70,7 @@ def initialize(self, Xq, Yq): """ q = len(Xq[0,:]) Xqhat, Yqhat = np.zeros(Xq.shape), np.zeros(Yq.shape) - if self.timestep == 0 and self.n <= q: + if self.timestep == 0 and np.linalg.matrix_rank(Xq) == self.n: weight = np.sqrt(self.weighting)**range(q-1,-1,-1) Xqhat, Yqhat = weight*Xq, weight*Yq self.A = Yqhat.dot(np.linalg.pinv(Xqhat)) diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index d889614..b9c179d 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -21,20 +21,22 @@ class WindowDMD: Usage: wdmd = WindowDMD(n,windowsize) - wdmd.initialize(Xq,Yq) + wdmd.initialize(Xw,Yw) wdmd.update(xold,yold,xnew,ynew) evals, modes = wdmd.computemodes() properties: n: state dimension windowsize: window size + weighting: weighting factor timestep: number of snapshot pairs processed (i.e., the current time step) A: DMD matrix, size n by n P: Matrix that contains information about recent w snapshots, size n by n methods: - initialize(Xq, Yq), initialize window DMD algorithm - update(x,y), update DMD computation when new snapshot pair (x,y) becomes available + initialize(Xw, Yw), initialize window DMD algorithm with w snapshot pairs + update(xold,yold,xnew,ynew), update DMD computation by forgetting old snapshot pair + and remember new snapshot pair computemodes(), compute and return DMD eigenvalues and DMD modes Authors: @@ -51,13 +53,14 @@ class WindowDMD: To import the WindowDMD class, add import window at head of Python scripts. To look up this documentation, type help(window.WindowDMD) or window.WindowDMD? """ - def __init__(self, n=0, windowsize=0, timestep=0, A=None, P=None): + def __init__(self, n=0, windowsize=0, weighting=1, timestep=0, A=None, P=None): """ Creat an object for window DMD - Usage: wdmd = WindowDMD(n,windowsize) + Usage: wdmd = WindowDMD(n,windowsize,weighting) """ self.n = n self.windowsize = windowsize + self.weighting = weighting self.timestep = timestep if A is None or P is None: self.A = np.zeros([n,n]) @@ -66,14 +69,16 @@ def __init__(self, n=0, windowsize=0, timestep=0, A=None, P=None): self.A = A self.P = P - def initialize(self, Xq, Yq): - """Initialize window DMD with first q snapshot pairs stored in (Xq, Yq) - Usage: wdmd.initialize(Xq,Yq) + def initialize(self, Xw, Yw): + """Initialize window DMD with first w snapshot pairs stored in (Xw, Yw) + Usage: wdmd.initialize(Xw,Yw) """ - q = len(Xq[0,:]) - if self.timestep == 0 and self.windowsize == q and self.windowsize >= self.n + 1: - self.A = Yq.dot(np.linalg.pinv(Xq)) - self.P = np.linalg.inv(Xq.dot(Xq.T)) + q = len(Xw[0,:]) + if self.timestep == 0 and self.windowsize == q and np.linalg.matrix_rank(Xw) == self.n: + weight = np.sqrt(self.weighting)**range(q-1,-1,-1) + Xwhat, Ywhat = weight*Xw, weight*Yw + self.A = Ywhat.dot(np.linalg.pinv(Xwhat)) + self.P = np.linalg.inv(Xwhat.dot(Xwhat.T))/self.weighting self.timestep += q def update(self, xold, yold, xnew, ynew): @@ -89,17 +94,17 @@ def update(self, xold, yold, xnew, ynew): # direct rank-2 update # define matrices U, V = np.column_stack((xold, xnew)), np.column_stack((yold, ynew)) - C = np.diag([-1,1]) + C = np.diag([-(self.weighting)**(self.windowsize),1]) # compute PkU matrix vector product beforehand PkU = self.P.dot(U) # compute AkU matrix vector product beforehand AkU = self.A.dot(U) # compute Gamma - Gamma = np.linalg.inv(C+U.T.dot(PkU)) + Gamma = np.linalg.inv(np.linalg.inv(C)+U.T.dot(PkU)) # update A self.A += (V-AkU).dot(Gamma).dot(PkU.T) # update P - self.P += -PkU.dot(Gamma).dot(PkU.T) + self.P = (self.P - PkU.dot(Gamma).dot(PkU.T))/self.weighting # ensure P is SPD by taking its symmetric part self.P = (self.P + self.P.T)/2 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 7812dd0..47a342d 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -104,13 +104,13 @@ def dyn(x,t): # Online DMD, weighting = 0.9 -evalsonlineDMD09 = np.empty((n,m),dtype=complex) +evalsonlineDMD2 = np.empty((n,m),dtype=complex) odmd = dmdtools.OnlineDMD(n,0.9) odmd.initialize(x[:,:q],y[:,:q]) start = time.clock() for k in range(q,m): odmd.update(x[:,k],y[:,k]) - evalsonlineDMD09[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt + evalsonlineDMD2[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt end = time.clock() print "Online DMD, weighting = 0.9, time = " + str(end-start) + " secs" @@ -122,7 +122,7 @@ def dyn(x,t): plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='Batch',linewidth=2.0) plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='Online, $wf=1$',linewidth=2.0) -plt.plot(t[q:], np.imag(evalsonlineDMD09[0,q:]), 'b-',label='Online, $wf=0.9$',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD2[0,q:]), 'b--',label='Online, $wf=0.9$',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 8b4dee4..b6ec452 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -90,19 +90,31 @@ def dyn(x,t): AminibatchDMD[:,:,k] = y[:,k-w+1:k+1].dot(np.linalg.pinv(x[:,k-w+1:k+1])) evalsminibatchDMD[:,k] = np.log(np.linalg.eigvals(AminibatchDMD[:,:,k]))/dt end = time.clock() -print "Mini-batch DMD, time = " + str(end-start) + " secs" +print "Mini-batch DMD, w = 10, time = " + str(end-start) + " secs" -# Window DMD, w = 10 -evalswindowDMD = np.empty((n,m),dtype=complex) -wdmd = dmdtools.WindowDMD(n,w) +# Window DMD, w = 10, weighting = 1 +evalswindowDMD1 = np.empty((n,m),dtype=complex) +wdmd = dmdtools.WindowDMD(n,w,1) wdmd.initialize(x[:,:w],y[:,:w]) start = time.clock() for k in range(w,m): wdmd.update(x[:,k-w],y[:,k-w],x[:,k],y[:,k]) - evalswindowDMD[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt + evalswindowDMD1[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt end = time.clock() -print "Window DMD, time = " + str(end-start) + " secs" +print "Window DMD, w=10, weighting = 1, time = " + str(end-start) + " secs" + + +# Window DMD, w = 10, weighting = 0.5 +evalswindowDMD2 = np.empty((n,m),dtype=complex) +wdmd = dmdtools.WindowDMD(n,w,0.5) +wdmd.initialize(x[:,:w],y[:,:w]) +start = time.clock() +for k in range(w,m): + wdmd.update(x[:,k-w],y[:,k-w],x[:,k],y[:,k]) + evalswindowDMD2[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt +end = time.clock() +print "Window DMD, w=10, weighting = 0.5, time = " + str(end-start) + " secs" # visualize true, batch, window @@ -111,7 +123,8 @@ def dyn(x,t): plt.rc('font', family='serif') plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, $w=10$',linewidth=2.0) -plt.plot(t[w:], np.imag(evalswindowDMD[0,w:]), 'g--',label='Window, $w=10$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD1[0,w:]), 'g--',label='Window, $w=10$, $wf=1$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD2[0,w:]), 'b--',label='Window, $w=10$, $wf=0.5$',linewidth=2.0) plt.tick_params(labelsize=20) plt.xlabel('Time', fontsize=20) plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) From b469ac7ba3df0d4ac1eaaba979e2e62bc2c237db Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Thu, 1 Jun 2017 13:45:15 -0400 Subject: [PATCH 27/32] update space complexity of windowed DMD --- matlab/OnlineDMD.m | 6 ++++-- matlab/WindowDMD.m | 8 +++++--- python/dmdtools/online.py | 5 ++++- python/dmdtools/window.py | 7 +++++-- 4 files changed, 18 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 59286e5..8088e88 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -1,5 +1,5 @@ % OnlineDMD is a class that implements online dynamic mode decomposition -% The time complexity (for one iteration) is O(n^2), and space complexity is +% The time complexity (multiply?add operation for one iteration) is O(n^2), and space complexity is % O(n^2), where n is the state dimension. % % Algorithm description: @@ -8,9 +8,11 @@ % 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) +% 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) diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 939b82b..c229355 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -1,6 +1,6 @@ % WindowDMD is a class that implements window dynamic mode decomposition -% The time complexity (for one iteration) is O(n^2), and space complexity is -% O(n^2), where n is the state dimension +% The time complexity (multiply?add operation for one iteration) is O(n^2), and space complexity is +% O(wn+2n^2), where n is the state dimension, w is the window size. % % Algorithm description: % At time step k, define two matrix @@ -14,6 +14,8 @@ % 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) @@ -75,7 +77,7 @@ end function initialize(obj, Xw, Yw) - % Initialize WnlineDMD with w snapshot pairs stored in (Xw, Yw) + % Initialize WindowDMD with w snapshot pairs stored in (Xw, Yw) % Usage: wdmd.initialize(Xw,Yw) q = length(Xw(1,:)); if(obj.timestep == 0 && obj.w == q && rank(Xw) == obj.n) diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 4dfb5bc..f3daf1d 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -1,9 +1,10 @@ +# -*- coding: utf-8 -*- import numpy as np class OnlineDMD: """OnlineDMD is a class that implements online dynamic mode decomposition - The time complexity (for one iteration) is O(n^2), and space complexity is + The time complexity (multiply–add operation for one iteration) is O(n^2), and space complexity is O(n^2), where n is the state dimension. Algorithm description: @@ -14,6 +15,8 @@ class OnlineDMD: should be measurements correponding to consecutive states z(k-1) and z(k). We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algrithm. + An exponential weighting factor can be used to place more weight on + recent data. Usage: odmd = OnlineDMD(n,weighting) diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index b9c179d..04709f5 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -1,10 +1,11 @@ +# -*- coding: utf-8 -*- import numpy as np class WindowDMD: """WindowDMD is a class that implements window dynamic mode decomposition - The time complexity (for one iteration) is O(n^2), and space complexity is - O(n^2), where n is the state dimension + The time complexity (multiply–add operation for one iteration) is O(n^2), and space complexity is + O(wn+2n^2), where n is the state dimension, w is the window size. Algorithm description: At time step k, define two matrix @@ -18,6 +19,8 @@ class WindowDMD: and remember the 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 algrithm. + An exponential weighting factor can be used to place more weight on + recent data. Usage: wdmd = WindowDMD(n,windowsize) From dab41f3cbb2cbf49bbb1d7d441a9b96d8b635342 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Thu, 1 Jun 2017 13:52:39 -0400 Subject: [PATCH 28/32] fix typo --- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 2 +- python/dmdtools/online.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 8088e88..31df361 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -1,5 +1,5 @@ % OnlineDMD is a class that implements online dynamic mode decomposition -% The time complexity (multiply?add operation for one iteration) is O(n^2), and space complexity is +% The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is % O(n^2), where n is the state dimension. % % Algorithm description: diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index c229355..eeeba3b 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -1,5 +1,5 @@ % WindowDMD is a class that implements window dynamic mode decomposition -% The time complexity (multiply?add operation for one iteration) is O(n^2), and space complexity is +% The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is % O(wn+2n^2), where n is the state dimension, w is the window size. % % Algorithm description: diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index f3daf1d..5635a9f 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -4,7 +4,7 @@ class OnlineDMD: """OnlineDMD is a class that implements online dynamic mode decomposition - The time complexity (multiply–add operation for one iteration) is O(n^2), and space complexity is + The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is O(n^2), where n is the state dimension. Algorithm description: From ca9c355825131f2343872d451b943c77809a9ada Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Wed, 21 Jun 2017 14:31:32 -0400 Subject: [PATCH 29/32] update window DMD implementation --- matlab/OnlineDMD.m | 38 +++++++++------- matlab/WindowDMD.m | 80 ++++++++++++++++++++++------------ matlab/online_demo.m | 36 +++++++-------- matlab/window_demo.m | 49 ++++++++++----------- python/dmdtools/online.py | 50 +++++++++++---------- python/dmdtools/window.py | 82 ++++++++++++++++++++--------------- python/scripts/online_demo.py | 27 ++++++------ python/scripts/window_demo.py | 32 +++++++------- 8 files changed, 219 insertions(+), 175 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 31df361..0959042 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -1,13 +1,15 @@ % OnlineDMD is a class that implements online dynamic mode decomposition -% The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is -% O(n^2), where n is the state dimension. +% 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, 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. @@ -23,19 +25,21 @@ % % properties: % n: state dimension -% weighting: weighting factor between 0 and 1 +% 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 +% 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 mdoes +% computemodes(), compute and return DMD eigenvalues and DMD modes % % Authors: % Hao Zhang @@ -87,7 +91,8 @@ function initialize(obj, Xq, Yq) end function initializeghost(obj) - % Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs before t=0 + % 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); @@ -96,11 +101,12 @@ function initializeghost(obj) 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). + % 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 + % compute P*x matrix vector product beforehand Px = obj.P*x; % Compute gamma gamma = 1/(1+x'*Px); @@ -115,7 +121,7 @@ function update(obj, x, y) end function [evals, modes] = computemodes(obj) - % Compute and return DMD eigenvalues and DMD modes at current time + % Compute DMD eigenvalues and DMD modes at current time % Usage: [evals, modes] = odmd.modes() [modes, evals] = eig(obj.A, 'vector'); end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index eeeba3b..d82bda1 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -1,17 +1,20 @@ % WindowDMD is a class that implements window dynamic mode decomposition -% The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is -% O(wn+2n^2), where n is the state dimension, w is the window size. +% 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, +% 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) +% 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 @@ -20,24 +23,28 @@ % Usage: % wdmd = WindowDMD(n,w,weighting) % wdmd.initialize(Xw,Yw) -% wdmd.update(xold, yold, xnew, ynew) +% wdmd.update(xnew, ynew) % [evals, modes] = wdmd.computemodes() % % properties: % n: state dimension % w: finite time window size -% weighting: weighting factor +% 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 +% P: Matrix that contains information about recent w snapshots, +% size n by n % % methods: % initialize(Xw, Yw), initialize window DMD algorithm -% update(xold, yold, xnew, ynew), update by forgetting old snapshot pairs, -% and remeber new snapshot pair -% At time k+1, if 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 xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) -% computemodes(), compute and return DMD eigenvalues and DMD mdoes +% 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 @@ -59,8 +66,11 @@ w = 0; % window size weighting = 1; % weighting factor 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 + P; % Matrix that contains information about recent w snapshots + % , size n by n end methods @@ -71,6 +81,8 @@ obj.n = n; obj.w = w; obj.weighting = weighting; + obj.Xw = zeros(n,w); + obj.Yw = zeros(n,w); obj.A = zeros(n,n); obj.P = zeros(n,n); end @@ -79,6 +91,10 @@ 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); @@ -90,16 +106,24 @@ function initialize(obj, Xw, Yw) obj.timestep = obj.timestep + q; end - function update(obj, xold, yold, xnew, ynew) + 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) - % Update the DMD computation by sliding the finite time window forward - % Forget the oldest pair of snapshots (xold, yold), and remebers 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 - % xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) - % Usage: wdmd.update(xold, yold, 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 @@ -123,7 +147,7 @@ function update(obj, xold, yold, xnew, ynew) end function [evals, modes] = computemodes(obj) - % Compute and return DMD eigenvalues and DMD modes at current time step + % Compute DMD eigenvalues and DMD modes at current time step % Usage: [evals, modes] = wdmd.computemodes() [modes, evals] = eig(obj.A, 'vector'); end diff --git a/matlab/online_demo.m b/matlab/online_demo.m index b2f03c0..8a25dc6 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -3,29 +3,29 @@ % We take a 2D time varying system given by dx/dt = A(t)x % where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], % w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) -% are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit +% are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. % -% 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, we would like to compute -% Ak = Yk*pinv(Xk). At time step K+1, we need to include new snapshot pair x(k+1), y(k+1) -% This can be done by brute-force batch DMD, -% and by efficient rank-1 updating online DMD algrithm. +% 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, +% we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force +% batch DMD, and by efficient rank-1 updating online DMD algrithm. Batch DMD +% computes DMD matrix by brute-force taking the pseudo-inverse directly. +% Online DMD computes the DMD matrix by using efficient rank-1 update idea. % -% Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly +% We compare the performance of online DMD (with weighting=1,0.9) with the +% brute-force batch DMD approach in terms of tracking time varying eigenvalues, +% by comparison with the analytical solution. Online DMD (weighting=1) and +% batch DMD should agree with each other (up to machine round-offer errors). % -% Online DMD computes the DMD matrix by using efficient rank-1 update idea -% We compare the performance of online DMD (with alpha=1,0.9) with the brute-force batch DMD -% approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -% Online DMD (weighting=1) and batch DMD should agree with each other (up to machine round-offer errors) % 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", -% in production, 2017. Available on arXiv. +% Hao Zhang +% Clarence W. Rowley % +% References: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Online Dynamic Mode Decomposition for Time-varying Systems", +% in production, 2017. Available on arXiv. +% % Date created: April 2017 % define dynamics diff --git a/matlab/window_demo.m b/matlab/window_demo.m index f9967f6..88f0b75 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -3,33 +3,32 @@ % We take a 2D time varying system given by dx/dt = A(t)x % where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], % w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) -% are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit +% are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. % -% 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, -% we would like to compute Ak = Yk*pinv(Xk). -% 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) -% This can be done by brute-force mini-batch DMD, -% and by efficient rank-2 updating window DMD algrithm. -% -% Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly -% -% Window DMD computes the DMD matrix by using efficient rank-2 update idea +% 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, we would like to compute Ak = Yk*pinv(Xk). This can +% be done by brute-force mini-batch DMD, and by efficient rank-2 updating window +% DMD algrithm. For window DMD, at time k+1, we need to forget the old snapshot +% pair xold = x(k-w+1), yold = y(k-w+1), and remember the new snapshot pair xnew +% = x(k+1), ynew = y(k+1). Mini-batch DMD computes DMD matrix by taking the +% pseudo-inverse directly. Window DMD computes the DMD matrix by using efficient +% rank-2 update idea. % % We compare the performance of window DMD with the brute-force mini-batch DMD -% approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -% They should agree with each other (up to machine round-offer errors) -% +% approach in terms of tracking time varying eigenvalues, by comparison with +% the analytical solution. They should agree with each other (up to machine +% round-offer errors). +% % 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", -% in production, 2017. Available on arXiv. -% +% Hao Zhang +% Clarence W. Rowley +% +% References: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% ``Online Dynamic Mode Decomposition for Time-varying Systems", +% in production, 2017. Available on arXiv. +% % Date created: April 2017 % define dynamics @@ -86,7 +85,7 @@ % window DMD tic for k = w+1:m - wdmd.update(x(:,k-w), y(:,k-w), x(:,k), y(:,k)); + wdmd.update(x(:,k), y(:,k)); evalswindowDMD1(:,k) = log(eig(wdmd.A))/dt; end elapsed_time = toc; @@ -100,7 +99,7 @@ % window DMD tic for k = w+1:m - wdmd.update(x(:,k-w), y(:,k-w), x(:,k), y(:,k)); + wdmd.update(x(:,k), y(:,k)); evalswindowDMD2(:,k) = log(eig(wdmd.A))/dt; end elapsed_time = toc; diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 5635a9f..d04e847 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -4,15 +4,17 @@ class OnlineDMD: """OnlineDMD is a class that implements online dynamic mode decomposition - The time complexity (multiply-add operation for one iteration) is O(n^2), and space complexity is - O(n^2), where n is the state dimension. + 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 correponding to consecutive states z(k-1) and z(k). + 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 correponding to consecutive + states z(k-1) and z(k). We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively by efficient rank-1 updating online DMD algrithm. An exponential weighting factor can be used to place more weight on @@ -28,14 +30,17 @@ class OnlineDMD: properties: n: state dimension weighting: weighting factor between (0,1] - timestep: number of snapshot pairs processed (i.e., the current time step) + timestep: number of snapshot pairs processed (i.e., current time step) 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 first 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 DMD computation when new snapshot pair (x,y) becomes available + initialize(Xq, Yq), initialize online DMD algorithm with first 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 DMD computation when new snapshot pair (x,y) + becomes available computemodes(), compute and return DMD eigenvalues and DMD modes Authors: @@ -50,22 +55,19 @@ class OnlineDMD: Date created: April 2017 To import the OnlineDMD class, add import online at head of Python scripts. - To look up this documentation, type help(online.OnlineDMD) or online.OnlineDMD? + To look up this documentation, type help(online.OnlineDMD) or + online.OnlineDMD? """ - def __init__(self, n=0, weighting=1, timestep=0, A=None, P=None): + def __init__(self, n=0, weighting=1): """ Creat an object for online DMD Usage: odmd = OnlineDMD(n,weighting) """ self.n = n self.weighting = weighting - self.timestep = timestep - if A is None or P is None: - self.A = np.zeros([n,n]) - self.P = np.zeros([n,n]) - else: - self.A = A - self.P = P + self.timestep = 0 + self.A = np.zeros([n,n]) + self.P = np.zeros([n,n]) def initialize(self, Xq, Yq): """Initialize online DMD with first q snapshot pairs stored in (Xq, Yq) @@ -81,7 +83,8 @@ def initialize(self, Xq, Yq): self.timestep += q def initializeghost(self): - """Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs before t=0 + """Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs + before t=0 Usage: odmd.initilizeghost() """ epsilon=1e-15 @@ -91,8 +94,9 @@ def initializeghost(self): def update(self, 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). + 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 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 04709f5..95a61aa 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -4,42 +4,45 @@ class WindowDMD: """WindowDMD is a class that implements window dynamic mode decomposition - The time complexity (multiply–add operation for one iteration) is O(n^2), and space complexity is - O(wn+2n^2), where n is the state dimension, w is the window size. + 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 correponding to consecutive states z(k-1) and z(k). - At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), - and remember the new snapshot pair xnew = x(k+1), ynew = y(k+1) + 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 correponding to consecutive + states z(k-1) and z(k). + At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), + yold = y(k-w+1), and remember the 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 algrithm. An exponential weighting factor can be used to place more weight on recent data. Usage: - wdmd = WindowDMD(n,windowsize) + wdmd = WindowDMD(n,w) wdmd.initialize(Xw,Yw) - wdmd.update(xold,yold,xnew,ynew) + wdmd.update(xnew,ynew) evals, modes = wdmd.computemodes() properties: n: state dimension - windowsize: window size + w: window size weighting: weighting factor - timestep: number of snapshot pairs processed (i.e., the current time step) + timestep: number of snapshot pairs processed (i.e., current time step) + 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, 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 with w snapshot pairs - update(xold,yold,xnew,ynew), update DMD computation by forgetting old snapshot pair - and remember new snapshot pair + update(xnew,ynew), update DMD computation by adding a new snapshot pair computemodes(), compute and return DMD eigenvalues and DMD modes Authors: @@ -54,50 +57,59 @@ class WindowDMD: Date created: April 2017 To import the WindowDMD class, add import window at head of Python scripts. - To look up this documentation, type help(window.WindowDMD) or window.WindowDMD? + To look up this documentation, type help(window.WindowDMD) + or window.WindowDMD? """ - def __init__(self, n=0, windowsize=0, weighting=1, timestep=0, A=None, P=None): + def __init__(self, n=0, w=0, weighting=1): """ Creat an object for window DMD - Usage: wdmd = WindowDMD(n,windowsize,weighting) + Usage: wdmd = WindowDMD(n,w,weighting) """ self.n = n - self.windowsize = windowsize + self.w = w self.weighting = weighting - self.timestep = timestep - if A is None or P is None: - self.A = np.zeros([n,n]) - self.P = np.zeros([n,n]) - else: - self.A = A - self.P = P + self.timestep = 0 + self.Xw = np.zeros([n,w]) + self.Yw = np.zeros([n,w]) + self.A = np.zeros([n,n]) + self.P = np.zeros([n,n]) def initialize(self, Xw, Yw): """Initialize window DMD with first w snapshot pairs stored in (Xw, Yw) Usage: wdmd.initialize(Xw,Yw) """ + # initialize Xw, Yw + self.Xw, self.Yw = Xw, Yw + # initialize A, P q = len(Xw[0,:]) - if self.timestep == 0 and self.windowsize == q and np.linalg.matrix_rank(Xw) == self.n: + if self.timestep == 0 and self.w == q \ + and np.linalg.matrix_rank(Xw) == self.n: weight = np.sqrt(self.weighting)**range(q-1,-1,-1) Xwhat, Ywhat = weight*Xw, weight*Yw self.A = Ywhat.dot(np.linalg.pinv(Xwhat)) self.P = np.linalg.inv(Xwhat.dot(Xwhat.T))/self.weighting self.timestep += q - def update(self, xold, yold, xnew, ynew): + def update(self, 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 - xold = x(k-w+1), yold = y(k-w+1), xnew = x(k+1), ynew = y(k+1) - Usage: wdmd.update(xold, yold, xnew, ynew) + 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, yold = self.Xw[:,0], self.Yw[:,0] + # Update recent w snapshots + self.Xw = np.column_stack((self.Xw[:,1:], xnew)) + self.Yw = np.column_stack((self.Yw[:,1:], ynew)) + # direct rank-2 update # define matrices U, V = np.column_stack((xold, xnew)), np.column_stack((yold, ynew)) - C = np.diag([-(self.weighting)**(self.windowsize),1]) + C = np.diag([-(self.weighting)**(self.w),1]) # compute PkU matrix vector product beforehand PkU = self.P.dot(U) # compute AkU matrix vector product beforehand diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 47a342d..1d28dd1 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -4,20 +4,19 @@ We take a 2D time varying system given by dx/dt = A(t)x where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) -are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit - -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, we would like to compute -Ak = Yk*pinv(Xk). This can be done by brute-force batch DMD, -and by efficient rank-1 updating online DMD algrithm. - -Batch DMD computes DMD matrix by brute-force taking the pseudo-inverse directly - -Online DMD computes the DMD matrix by using efficient rank-1 update idea - -We compare the performance of online DMD (with weighting=1,0.9) with the brute-force batch DMD -approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -Online DMD (weighting=1) and batch DMD should agree with each other (up to machine round-offer errors) +are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. + +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, +we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force +batch DMD, and by efficient rank-1 updating online DMD algrithm. Batch DMD +computes DMD matrix by brute-force taking the pseudo-inverse directly. +Online DMD computes the DMD matrix by using efficient rank-1 update idea. + +We compare the performance of online DMD (with weighting=1,0.9) with the +brute-force batch DMD approach in terms of tracking time varying eigenvalues, +by comparison with the analytical solution. Online DMD (weighting=1) and +batch DMD should agree with each other (up to machine round-offer errors). Authors: Hao Zhang diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index b6ec452..74382bb 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -4,22 +4,22 @@ We take a 2D time varying system given by dx/dt = A(t)x where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) -are pure imaginary, i.e, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit - -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, -we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force mini-batch DMD, -and by efficient rank-2 updating window DMD algrithm. -For window DMD, at time k+1, we need to forget the old snapshot pair xold = x(k-w+1), yold = y(k-w+1), -and remember the new snapshot pair xnew = x(k+1), ynew = y(k+1) - -Mini-batch DMD computes DMD matrix by taking the pseudo-inverse directly - -Window DMD computes the DMD matrix by using efficient rank-2 update idea +are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. + +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, we would like to compute Ak = Yk*pinv(Xk). This can +be done by brute-force mini-batch DMD, and by efficient rank-2 updating window +DMD algrithm. For window DMD, at time k+1, we need to forget the old snapshot +pair xold = x(k-w+1), yold = y(k-w+1), and remember the new snapshot pair xnew += x(k+1), ynew = y(k+1). Mini-batch DMD computes DMD matrix by taking the +pseudo-inverse directly. Window DMD computes the DMD matrix by using efficient +rank-2 update idea. We compare the performance of window DMD with the brute-force mini-batch DMD -approach in terms of tracking time varying eigenvalues, by comparison with the analytical solution -They should agree with each other (up to machine round-offer errors) +approach in terms of tracking time varying eigenvalues, by comparison with +the analytical solution. They should agree with each other (up to machine +round-offer errors). Authors: Hao Zhang @@ -99,7 +99,7 @@ def dyn(x,t): wdmd.initialize(x[:,:w],y[:,:w]) start = time.clock() for k in range(w,m): - wdmd.update(x[:,k-w],y[:,k-w],x[:,k],y[:,k]) + wdmd.update(x[:,k],y[:,k]) evalswindowDMD1[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt end = time.clock() print "Window DMD, w=10, weighting = 1, time = " + str(end-start) + " secs" @@ -111,7 +111,7 @@ def dyn(x,t): wdmd.initialize(x[:,:w],y[:,:w]) start = time.clock() for k in range(w,m): - wdmd.update(x[:,k-w],y[:,k-w],x[:,k],y[:,k]) + wdmd.update(x[:,k],y[:,k]) evalswindowDMD2[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt end = time.clock() print "Window DMD, w=10, weighting = 0.5, time = " + str(end-start) + " secs" From a580e1e80ca9890998b8decaf892bcb6479ecdf9 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Wed, 21 Jun 2017 17:05:48 -0400 Subject: [PATCH 30/32] minor change to documentation --- matlab/OnlineDMD.m | 3 ++- matlab/WindowDMD.m | 7 ++++--- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 6 +++--- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 0959042..7986ba9 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -58,7 +58,7 @@ classdef OnlineDMD < handle properties n = 0; % state dimension - weighting = 1; % weighting factor + weighting = 1; % weighting factor in (0,1] timestep = 0; % number of snapshots processed A; % DMD matrix P; % matrix that contains information about past snapshots @@ -71,6 +71,7 @@ if nargin == 2 obj.n = n; obj.weighting = weighting; + obj.timestep = 0; obj.A = zeros(n,n); obj.P = zeros(n,n); end diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index d82bda1..0f4cc1a 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -64,7 +64,7 @@ properties n = 0; % state dimension w = 0; % window size - weighting = 1; % weighting factor + 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 @@ -81,6 +81,7 @@ 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); @@ -129,9 +130,9 @@ function update(obj, xnew, ynew) % define matrices U = [xold, xnew]; V = [yold, ynew]; C = diag([-(obj.weighting)^(obj.w),1]); - % compute PkU matrix vector product beforehand + % compute PkU matrix matrix product beforehand PkU = obj.P*U; - % compute AkU matrix vector product beforehand + % compute AkU matrix matrix product beforehand AkU = obj.A*U; % compute Gamma Gamma = inv(inv(C)+U'*PkU); diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index d04e847..cf18c80 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -29,7 +29,7 @@ class OnlineDMD: properties: n: state dimension - weighting: weighting factor between (0,1] + weighting: weighting factor in (0,1] timestep: number of snapshot pairs processed (i.e., current time step) A: DMD matrix, size n by n P: Matrix that contains information about past snapshots, size n by n diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 95a61aa..266511e 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -33,7 +33,7 @@ class WindowDMD: properties: n: state dimension w: window size - weighting: weighting factor + weighting: weighting factor in (0,1] timestep: number of snapshot pairs processed (i.e., current time step) Xw: recent w snapshots x stored in Xw, size n by w Yw: recent w snapshots y stored in Yw, size n by w @@ -110,9 +110,9 @@ def update(self, xnew, ynew): # define matrices U, V = np.column_stack((xold, xnew)), np.column_stack((yold, ynew)) C = np.diag([-(self.weighting)**(self.w),1]) - # compute PkU matrix vector product beforehand + # compute PkU matrix matrix product beforehand PkU = self.P.dot(U) - # compute AkU matrix vector product beforehand + # compute AkU matrix matrix product beforehand AkU = self.A.dot(U) # compute Gamma Gamma = np.linalg.inv(np.linalg.inv(C)+U.T.dot(PkU)) From 6fe3f2defe9cd1fae2011c5f363ee1d96bd02d36 Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Mon, 10 Jul 2017 21:09:52 -0400 Subject: [PATCH 31/32] update reference to include arXiv number --- matlab/OnlineDMD.m | 4 ++-- matlab/WindowDMD.m | 4 ++-- matlab/online_demo.m | 6 +++--- matlab/window_demo.m | 6 +++--- python/dmdtools/online.py | 4 ++-- python/dmdtools/window.py | 4 ++-- python/scripts/online_demo.py | 5 +++-- python/scripts/window_demo.py | 5 +++-- 8 files changed, 20 insertions(+), 18 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index 7986ba9..c36f618 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -47,8 +47,8 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. Available on arXiv. +% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% arXiv preprint arXiv:1707.02876, 2017. % % Created: % April 2017. diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 0f4cc1a..256fd25 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -52,8 +52,8 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. Available on arXiv. +% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% arXiv preprint arXiv:1707.02876, 2017. % % Created: % April 2017. diff --git a/matlab/online_demo.m b/matlab/online_demo.m index 8a25dc6..b849dfc 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -22,9 +22,9 @@ % Clarence W. Rowley % % References: -% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. Available on arXiv. +% 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. % % Date created: April 2017 diff --git a/matlab/window_demo.m b/matlab/window_demo.m index 88f0b75..f25cd42 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -25,9 +25,9 @@ % Clarence W. Rowley % % References: -% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems", -% in production, 2017. Available on arXiv. +% 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. % % Date created: April 2017 diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index cf18c80..2ca503a 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -49,8 +49,8 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. Available on arXiv. + ``Online Dynamic Mode Decomposition for Time-varying Systems,” + arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 266511e..90a40ee 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -51,8 +51,8 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. Available on arXiv. + ``Online Dynamic Mode Decomposition for Time-varying Systems,” + arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 1d28dd1..2fe6f34 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """ An example to demonstrate online dynamic mode decomposition @@ -24,8 +25,8 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. Available on arXiv. + ``Online Dynamic Mode Decomposition for Time-varying Systems,” + arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 """ diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index 74382bb..aaf77e3 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """ An example to demonstrate window dynamic mode decomposition @@ -27,8 +28,8 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems", - in production, 2017. Available on arXiv. + ``Online Dynamic Mode Decomposition for Time-varying Systems,” + arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 """ From 3a663af55a83136328e4684fdcfba244777924bd Mon Sep 17 00:00:00 2001 From: Hao Zhang Date: Mon, 10 Jul 2017 21:27:45 -0400 Subject: [PATCH 32/32] update references --- matlab/OnlineDMD.m | 2 +- matlab/WindowDMD.m | 2 +- matlab/online_demo.m | 2 +- matlab/window_demo.m | 2 +- python/dmdtools/online.py | 2 +- python/dmdtools/window.py | 2 +- python/scripts/online_demo.py | 2 +- python/scripts/window_demo.py | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m index c36f618..9d82f6e 100644 --- a/matlab/OnlineDMD.m +++ b/matlab/OnlineDMD.m @@ -47,7 +47,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% "Online Dynamic Mode Decomposition for Time-varying Systems," % arXiv preprint arXiv:1707.02876, 2017. % % Created: diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m index 256fd25..75b95a3 100644 --- a/matlab/WindowDMD.m +++ b/matlab/WindowDMD.m @@ -52,7 +52,7 @@ % % Reference: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% "Online Dynamic Mode Decomposition for Time-varying Systems," % arXiv preprint arXiv:1707.02876, 2017. % % Created: diff --git a/matlab/online_demo.m b/matlab/online_demo.m index b849dfc..d2fa825 100644 --- a/matlab/online_demo.m +++ b/matlab/online_demo.m @@ -23,7 +23,7 @@ % % References: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% "Online Dynamic Mode Decomposition for Time-varying Systems," % arXiv preprint arXiv:1707.02876, 2017. % % Date created: April 2017 diff --git a/matlab/window_demo.m b/matlab/window_demo.m index f25cd42..2f81830 100644 --- a/matlab/window_demo.m +++ b/matlab/window_demo.m @@ -26,7 +26,7 @@ % % References: % Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, -% ``Online Dynamic Mode Decomposition for Time-varying Systems,? +% "Online Dynamic Mode Decomposition for Time-varying Systems," % arXiv preprint arXiv:1707.02876, 2017. % % Date created: April 2017 diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py index 2ca503a..fd79138 100644 --- a/python/dmdtools/online.py +++ b/python/dmdtools/online.py @@ -49,7 +49,7 @@ class OnlineDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems,” + "Online Dynamic Mode Decomposition for Time-varying Systems," arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py index 90a40ee..2d072fd 100644 --- a/python/dmdtools/window.py +++ b/python/dmdtools/window.py @@ -51,7 +51,7 @@ class WindowDMD: References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems,” + "Online Dynamic Mode Decomposition for Time-varying Systems," arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py index 2fe6f34..6a824c4 100644 --- a/python/scripts/online_demo.py +++ b/python/scripts/online_demo.py @@ -25,7 +25,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems,” + "Online Dynamic Mode Decomposition for Time-varying Systems," arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017 diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py index aaf77e3..3ea0b20 100644 --- a/python/scripts/window_demo.py +++ b/python/scripts/window_demo.py @@ -28,7 +28,7 @@ References: Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, - ``Online Dynamic Mode Decomposition for Time-varying Systems,” + "Online Dynamic Mode Decomposition for Time-varying Systems," arXiv preprint arXiv:1707.02876, 2017. Date created: April 2017