# CCAMA – A customized alternating minimization algorithm for structured covariance completion problems

 Written by: Armin Zare and Mihailo R. Jovanovic April 2016 Matlab Files Download all Matlab files Papers Theory and algorithms:   arXiv:1412.3399 Application to turbulent flows:   arXiv:1602.05105

## Purpose

This website provides a Matlab implementation of a customized alternating minimization algorithm (CCAMA) for solving the covariance completion problem (CC). This problem aims to explain and complete partially available statistical signatures of dynamical systems by identifying the directionality and dynamics of input excitations into linear approximations of the system. In this, we seek to explain correlation data with the least number of possible input disturbance channels. This inverse problem is formulated as a rank minimization, and for its solution, we employ a convex relaxation based on the nuclear norm. CCAMA exploits the structure of (CC) in order to gain computational efficiency and improve scalability.

## Problem formulation

Consider a linear time-invariant system

where is the state vector, is the output, is a stationary zero-mean stochastic process, is the dynamic matrix, is the input matrix, and is the output matrix. For Hurwitz and controllable , a positive definite matrix qualifies as the steady-state covariance matrix of the state vector

if and only if the linear equation

is solvable for . Here, is the expectation operator, represents the cross-correlation of the state and the input , and denotes the complex conjugate transpose.

The algebraic relation between second-order statistics of the state and forcing can be used to explain partially known sampled second-order statistics using stochastically-driven LTI systems. While the dynamical generator is known, the origin and directionality of stochastic excitation is unknown. It is also important to restrict the complexity of the forcing model. This complexity is quantified by the number of degrees of freedom that are directly influenced by stochastic forcing and translates into the number of input channels or . It can be shown that the rank of is closely related to the signature of the matrix

The signature of a matrix is determined by the number of its positive, negative, and zero eigenvalues. In addition, the rank of bounds the rank of .

Based on this, the problem of identifying low-complexity structures for stochastic forcing can be formulated as the following structured covariance completion problem

Here, is a positive regularization parameter, the matrices , , , and are problem data, and the Hermitian matrices , are optimization variables. Entries of represent partially available second-order statistics of the output , the symbol denotes elementwise matrix multiplication, and is the structural identity matrix,

Convex optimization problem (CC) combines the nuclear norm with an entropy function in order to target low-complexity structures for stochastic forcing and facilitate construction of a particular class of low-pass filters that generate colored-in-time forcing correlations. The nuclear norm, i.e., the sum of singular values of a matrix, , is used as a proxy for rank minimization. On the other hand, the logarithmic barrier function in the objective is introduced to guarantee the positive definiteness of the state covariance matrix .

In (CC), determines the importance of the nuclear norm relative to the logarithmic barrier function. The convexity of (CC) follows from the convexity of the objective function

and the convexity of the constraint set. Problem (CC) can be equivalently expressed as follows,

where the constraints are now given by

Here, and are linear operators, with

and their adjoints, with respect to the standard inner product , are given by

By splitting into positive and negative definite parts,

it can be shown that (CC-) can be cast as an SDP,

The Lagrange dual of (P) is given by

where Hermitian matrices , are the dual variables associated with the equality constraints in (P) and is the number of states.

## Alternating Minimization Algorithm (AMA)

The logarithmic barrier function in (CC) is strongly convex over any compact subset of the positive definite cone. This makes it well-suited for the application of AMA, which requires strong convexity of the smooth part of the objective function.

The augmented Lagrangian associated with (CC-) is

where is a positive scalar and is the Frobenius norm.

AMA consists of the following steps:

These terminate when the duality gap

and the primal residual

are sufficiently small, i.e., and . In the -minimization step, AMA minimizes the Lagrangian with respect to . This step is followed by a -minimization step in which the augmented Lagrangian is minimized with respect to . Finally, the Lagrange multipliers, and , are updated based on the primal residuals with the step-size .

The -minimization step amounts to a matrix inversion and the -minimization step amounts to application of the soft-thresholding operator on the singular values of a matrix:

For Hermitian matrix with singular value decomposition , and denote the singular value thresholding and saturation operators, respectively:

with . The updates of Lagrange multipliers guarantee dual feasibility at each iteration, i.e., for all .

## Alternating Direction Method of Multipliers (ADMM)

In contrast to AMA, ADMM minimizes the augmented Lagrangian in each step of the iterative procedure. In addition, ADMM does not have efficient step-size selection rules. Typically, either a constant step-size is selected or the step-size is adjusted to keep the norms of primal and dual residuals within a constant factor of one another Boyd, et al. ’11.

ADMM consists of the following steps:

which terminate with similar stopping criteria to ADMM.

While the -minimization step is equivalent to that of AMA, the -update in ADMM is obtained by minimizing the augmented Lagrangian. This amounts to solving the following optimization problem

where and . From first order optimality conditions we have

Since and are not unitary operators, the -minimization step does not have an explicit solution. To solve this problem we use a proximal gradient method Parikh and Boyd ’13 to update . By linearizing the quadratic term around the current inner iterate and adding a quadratic penalty on the difference between and , is obtained as the minimizer of

To ensure convergence of the proximal gradient method, the parameter has to satisfy , where we use power iteration to compute the largest eigenvalue of the operator .

By taking the variation with respect to we arrive at the first order optimality condition

The solution to () can be iteratively computed as

where the th entry of the vector is given by

Here, 's are the eigenvalues of the matrix on the right-hand-side of () and is the matrix of the corresponding eigenvectors. As it is typically done in proximal gradient algorithms Parikh and Boyd ’13, starting with , we obtain by repeating inner iterations until the desired accuracy is reached.

## Acknowledgements

This project is supported by the