Description of dmdsp.m

% Sparsity-Promoting Dynamic Mode Decomposition (DMDSP)
%
% Written by Mihailo R. Jovanovic, August 2012
%
% The alternating direction method of multipliers (ADMM)
% is used to solve the sparsity-promoting DMD problem:
%
% minimize  J(x) + gamma * || x ||_1              (SP)
%
% J(x) = || G - L diag{x} R ||_F^2
%
% complex-valued data G, L, R; optimization variable x
%
% J(x) can be equivalently written as
% J(x) = x'*P*x - 2*real(q'*x) + s
%
% P = (L'*L).*conjugate(R*R')
% q = conjugate(diag(R*G'*L))
% s = trace(G'*G)
%
% For the identified sparsity pattern, the optimal vector of
% amplitudes is found as the answer to the following
% structured quadratic optimization problem:
%
%       minimize    J(x)
%                                                  (POL)
%       subject to  E'*x = 0
%
% where the columns of the matrix E are the unit vectors
% whose non-zero elements correspond to zero components of x
%
% Syntax:
%
% answer = dmdsp(P,q,s,gammaval,options)
%
% Inputs:  (1) matrix P
%              vector q
%              scalar s
%              sparsity promoting parameter gamma
%
%          (2) options
%
%              options.rho     - augmented Lagrangian parameter rho
%              options.maxiter - maximum number of ADMM iterations
%              options.eps_abs - absolute tolerance
%              options.eps_rel - relative tolerance
%
%              If options argument is omitted, the default values are set to
%
%              options.rho = 1
%              options.maxiter = 10000
%              options.eps_abs = 1.e-6
%              options.eps_rel = 1.e-4
%
% Output:  answer - gamma-parameterized structure containing
%
%          answer.gamma - sparsity-promoting parameter gamma
%          answer.xsp   - vector of amplitudes resulting from (SP)
%          answer.xpol  - vector of amplitudes resulting from (POL)
%          answer.Jsp   - J resulting from (SP)
%          answer.Jpol  - J resulting from (POL)
%          answer.Nz    - number of nonzero elements of x
%          answer.Ploss - optimal performance loss 100*sqrt(J(xpol)/J(0))
%
% Additional information:
%
% http://www.umn.edu/~mihailo/software/dmdsp/

function answer = dmdsp(P,q,s,gammaval,options)

% Initialization
if nargin < 4
    error('At least four input arguments are required.')
elseif nargin == 4
    options = struct('rho',1,'maxiter',10000,'eps_abs',1.e-6,'eps_rel',1.e-4);
elseif nargin > 5
    error('Too many input arguments.')
end

% Data preprocessing
rho = options.rho;
Max_ADMM_Iter = options.maxiter;
eps_abs = options.eps_abs;
eps_rel = options.eps_rel;

% Number of optimization variables
n = length(q);
% Identity matrix
I = eye(n);

% Allocate memory for gamma-dependent output variables
answer.gamma = gammaval;
answer.Nz    = zeros(1,length(gammaval)); % number of non-zero amplitudes
answer.Jsp   = zeros(1,length(gammaval)); % square of Frobenius norm (before polishing)
answer.Jpol  = zeros(1,length(gammaval)); % square of Frobenius norm (after polishing)
answer.Ploss = zeros(1,length(gammaval)); % optimal performance loss (after polishing)
answer.xsp   = zeros(n,length(gammaval)); % vector of amplitudes (before polishing)
answer.xpol  = zeros(n,length(gammaval)); % vector of amplitudes (after polishing)

% Cholesky factorization of matrix P + (rho/2)*I
Prho = (P + (rho/2)*I);
Plow = chol(Prho,'lower');
Plow_star = Plow';

for i = 1:length(gammaval),

    gamma = gammaval(i);

    % Initial conditions
    y = zeros(n,1); % Lagrange multiplier
    z = zeros(n,1); % copy of x

    % Use ADMM to solve the gamma-parameterized problem
    for ADMMstep = 1 : Max_ADMM_Iter,

        % x-minimization step
        u = z - (1/rho)*y;
        % x = (P + (rho/2)*I)\(q + (rho)*u)
        xnew = Plow_star\(Plow\(q + (rho/2)*u));

        % z-minimization step
        a = (gamma/rho)*ones(n,1);
        v = xnew + (1/rho)*y;
        % Soft-thresholding of v
        znew = ( (1 - a ./ abs(v)) .* v ) .* (abs(v) > a);

        % Primal and dual residuals
        res_prim = norm(xnew - znew);
        res_dual = rho * norm(znew - z);

        % Lagrange multiplier update step
        y = y + rho*(xnew - znew);

        % Stopping criteria
        eps_prim = sqrt(n) * eps_abs + eps_rel * max([norm(xnew),norm(znew)]);
        eps_dual = sqrt(n) * eps_abs + eps_rel * norm(y);

        if (res_prim < eps_prim) && (res_dual < eps_dual)
            break;
        else
            z = znew;
        end

    end

    % Record output data
    answer.xsp(:,i) = z; % vector of amplitudes
    answer.Nz(i) = nnz(answer.xsp(:,i)); % number of non-zero amplitudes
    answer.Jsp(i) = real(z'*P*z) - 2*real(q'*z) + s; % Frobenius norm (before polishing)

    % Polishing of the nonzero amplitudes
    % Form the constraint matrix E for E^T x = 0
    ind_zero = find( abs(z) < 1.e-12); % find indices of zero elements of z
    m = length(ind_zero); % number of zero elements
    E = I(:,ind_zero);
    E = sparse(E);

    % Form KKT system for the optimality conditions
    KKT = [P, E; E', zeros(m,m)];
    rhs = [q; zeros(m,1)];

    % Solve KKT system
    sol = KKT \ rhs;

    % Vector of polished (optimal) amplitudes
    xpol = sol(1:n);

    % Record output data
    answer.xpol(:,i) = xpol;
    % Polished (optimal) least-squares residual
    answer.Jpol(i) = real(xpol'*P*xpol) - 2*real(q'*xpol) + s;
    % Polished (optimal) performance loss
    answer.Ploss(i) = 100*sqrt(answer.Jpol(i)/s);

end

Published with MATLAB® 7.13