add stableRegression() member function

dmd.m

@@ -62,6 +62,52 @@ classdef dmd < handle
             % image space)
             obj.Phi = obj.X2*obj.V/obj.S*obj.W;
             
+            % tell user how long it took
+            regressionTime = toc;
+            if(obj.verbose)
+                fprintf('DMD regression() took %f seconds.\n',regressionTime);
+            end
+		end
+		
+		function stableRegression(obj,varargin) % argument is truncation rank r
+			% Follows Ian Manchester's not-yet published work on stable
+			% least squares
+            if nargin == 2
+                obj.r = varargin{1};
+            end
+            tic;
+            [obj.U0,obj.S0,obj.V0] = svd(obj.X,'econ');
+            
+            % truncate
+			if obj.r <= size(obj.S0,1)
+				obj.U = obj.U0(:,1:obj.r);
+				obj.S = obj.S0(1:obj.r,1:obj.r);
+				obj.V = obj.V0(:,1:obj.r);
+				fprintf('DMD: Truncated to rank %d out of %d.\n',obj.r,size(obj.S0,1));
+			else
+				fprintf('DMD: skipping truncation. Requested rank %d is higher than max possible rank %d.\n',obj.r,size(obj.S0,1));
+				obj.U = obj.U0;
+				obj.S = obj.S0;
+				obj.V = obj.V0;
+			end
+            
+            % construct weight-space dynamics
+			N = size(obj.X,2); % number of timesteps
+			xtilde = obj.U'*obj.X; % project into DMD "latent" state
+			ztilde = obj.U'*obj.X2; 
+			Sigma_xx = xtilde*xtilde'/N;
+			Sigma_zx = ztilde*xtilde'/N;
+			Sigma_zz = ztilde*ztilde'/N;
+			Sigma_s = (Sigma_xx + Sigma_zz)/2;
+			
+		
+			obj.Atilde = Sigma_zx/Sigma_s;
+            [obj.W,obj.Lambda] = eig(obj.Atilde);
+            
+            % Phi projects from weight-space to original data space (e.g.
+            % image space)
+            obj.Phi = obj.X2*obj.V/obj.S*obj.W;
+            
             % tell user how long it took
             regressionTime = toc;
             if(obj.verbose)