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

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

where

is the state vector, y(t)

is the output, u(t)

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, B)

, a positive definite matrix

qualifies as the steady-state covariance matrix of the state vector

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 {rm rank}(B)

. 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

$begin{array}{cl} {rm minimize} & -{rm log,det}left(Xright) ,+, gamma,||Z||_* [.25cm] {rm subject~to} & ~A , X ,+, X A^* ,+, Z ;=; 0 [.15cm] & ,,left(C X C^* right)circ E ,-, G ;=; 0. end{array} hspace{1.5cm} {rm (CC)}$

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 circ

denotes elementwise matrix multiplication, and

is the structural identity matrix,

$E_{ij} ;=; left{ begin{array}{ll} 1, &~ mathrm{if} ~ G_{ij} ~ {rm is~available} [.1cm] 0, &~ mathrm{if}~ G_{ij}~ {rm is~unavailable.} end{array} right.$

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, $||Z||_* mathrel{mathop:}= sum_{i} sigma_i (Z)$ , 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,

$begin{array}{cl} {rm minimize} & -{rm log,det} X ,+, gamma,||Z||_* [.25cm] {rm subject~to} & ~{cal A}, X ,+, {cal B}, Z ,-, {cal C} ,=, 0, end{array} hspace{1.5cm} {rm (CC}-{rm 1)}$

$left[ begin{array}{c} {cal A}_1 {cal A}_2 end{array} right] X ,+, left[ begin{array}{c} I 0 end{array} right] Z ,-, left[ begin{array}{c} 0 G end{array} right] ,=, 0.$

$begin{array}{c} {cal A}_1(X) ; mathrel{mathop:}= ; A , X ; + ; XA^* [.15cm] {cal A}_2(X) ; mathrel{mathop:}= ; (C , X , C^*)circ E end{array}$

and their adjoints, with respect to the standard inner product $left<M_1, M_2right> mathrel{mathop:}= {rm trace} , (M_1^* M_2)$ , are given by

$begin{array}{c} {cal A}_1^dagger (Y) ; = ; A^* , Y ;+; Y A [.15cm] {cal A}_2^dagger (Y) ; = ; C^* (Ecirc Y) , C. end{array}$

$begin{array}{cl} {rm minimize} & -{rm log,det} X ; + ; gamma left( {cal A}, (Z_+) , + , {rm trace}, (Z_-)right) [.25cm] {rm subject~to} & {cal A}_1(X) ,+, Z_+ , - , Z_- ,=, 0 [0.15cm] & {cal A}_2 (X) ,-, G ,=, 0 [0.15cm] & Z_+ , succeq , 0,~~~~Z_- , succeq , 0. end{array} hspace{1.5cm} {rm (P)}$

$begin{array}{cl} {rm maximize} & {rm log,det} left( {cal A}_1^dagger(Y_1) , + , {cal A}_2^dagger(Y_2) right) , - ; left<G,, Y_2 right> , + ; n [.25cm] {rm subject~to} & ||Y_1||_2 , leq , gamma end{array} hspace{1.5cm} {rm (D)}$

where Hermitian matrices Y_1

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.

$begin{array}{l} {cal L}_rho (X, Z; Y_1, Y_2) ; = ; displaystyle{-{rm log,det} X , + , gamma, ||Z||_*} , +, left<Y_1,, {cal A}_1 (X) + Z right> , + , left<Y_2,, {cal A}_2 (X) - G right> ~ + [0.25cm] hfill { displaystyle{ frac{rho}{2}, ||{cal A}_1 (X) , + , Z||_F^2} ; + ; displaystyle{ frac{rho}{2}, ||{cal A}_2 (X) , - , G||_F^2} } end{array}$

$begin{array}{rrcl} X-mathbf{minimization}~mathbf{step:} & ~~ X^{k+1} &!!mathrel{mathop:}=!!& {rm argmin}_{X} , {cal L}_0, ( X,, Z^k; , Y_1^k,, Y_2^k) [0.25cm] Z-mathbf{minimization}~mathbf{step:} & ~~ Z^{k+1} &!!mathrel{mathop:}=!!& {rm argmin}_{Z} , {cal L}_{rho}, ( X^{k+1},, Z; , Y_1^k,, Y_2^k) [0.25cm] Y-mathbf{update}~mathbf{step:} & ~~ Y_1^{k+1} &!!mathrel{mathop:}=!!& Y_1^k ,+, displaystyle{rho left( {cal A}_1 (X^{k+1}) , + , Z^{k+1} right)} [0.25cm] & ~~ Y_2^{k+1} &!!mathrel{mathop:}=!!& Y_2^k ,+, displaystyle{rho left( {cal A}_2 (X^{k+1}) , - , G right)} end{array}$

$Delta_{mathrm{gap}} ; mathrel{mathop:}= ; -{rm log,det} X^{k+1} ,+, gamma,||Z^{k+1}||_* , - , J_d left(Y_1^{k+1}, Y_2^{k+1} right)$

are sufficiently small, i.e., $|Delta_{mathrm{gap}}| leq epsilon_1$ and $Delta_{mathrm{p}} leq epsilon_2$ . In the

-minimization step, AMA minimizes the Lagrangian ${cal L}_0$ with respect to

. This step is followed by a

-minimization step in which the augmented Lagrangian ${cal L}_{rho}$ is minimized with respect to

. Finally, the Lagrange multipliers, Y_1

and

, are updated based on the primal residuals with the step-size rho

The

-minimization step amounts to a matrix inversion and the

-minimization step amounts to application of the soft-thresholding operator ${cal S}_tau$ on the singular values of a matrix:

$begin{array}{rcl} X^{k+1} &!!=!!& left( {cal A}_1^dagger (Y_1^k) ,+, {cal A}_2^dagger (Y_2^k) right)^{-1} [0.25cm] Z^{k+1} &!!=!!& {cal S}_tau left( - {cal A}_1 ( X^{k+1} ) , - , (1/rho)Y_1^k right), ~~~~ tau = gamma/rho [0.25cm] Y_1^{k+1} &!!=!!& {cal T}_{gamma} left(Y_1^k ,+, rho, {cal A}_1 (X^{k+1}) right) [0.25cm] Y_2^{k+1} &!!=!!& Y_2^k ,+, displaystyle{rho left( {cal A}_2 (X^{k+1}) , - , G right)} end{array}$

For Hermitian matrix

with singular value decomposition M = U Sigma U^*

, ${cal S}_tau$ and ${cal T}_tau$ denote the singular value thresholding and saturation operators, respectively:

$begin{array}{rclrcl} {cal S}_tau (M) &!! mathrel{mathop:}= !!& U , {cal S}_tau (Sigma) ,U^*, & ~~~~ {cal S}_tau (Sigma) &!! = !!& {rm diag} left( left(sigma_i ,-, tau right)_+ right) [.25cm] {cal T}_tau (M) &!! mathrel{mathop:}= !!& U , {cal T}_tau (Sigma) ,U^*, & ~~~~ {cal T}_tau (Sigma) & !! = !! & {rm diag} left( min left( max( sigma_i, -, tau ), tau right) right) end{array}$

with $a_+ mathrel{mathop:}= {rm max} , {a, 0}$ . The updates of Lagrange multipliers guarantee dual feasibility at each iteration, i.e., $||Y_1^{k+1}||_2 leq gamma$ 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.

$begin{array}{rrcl} X-mathbf{minimization}~mathbf{step:} & ~~ X^{k+1} &!!mathrel{mathop:}=!!& {rm argmin}_{X} , {cal L}_{rho}, ( X,, Z^k; , Y_1^k,, Y_2^k) [0.25cm] Z-mathbf{minimization}~mathbf{step:} & ~~ Z^{k+1} &!!mathrel{mathop:}=!!& {rm argmin}_{Z} , {cal L}_{rho}, ( X^{k+1},, Z; , Y_1^k,, Y_2^k) [0.25cm] Y-mathbf{update}~mathbf{step:} & ~~ Y_1^{k+1} &!!mathrel{mathop:}=!!& Y_1^k ,+, displaystyle{rho left( {cal A}_1 (X^{k+1}) , + , Z^{k+1} right)} [0.25cm] & ~~ Y_2^{k+1} &!!mathrel{mathop:}=!!& Y_2^k ,+, displaystyle{rho left( {cal A}_2 (X^{k+1}) , - , G right)} end{array}$

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

$begin{array}{cl} {rm minimize} &!! - {rm log,det} X ;+; displaystyle{frac{rho}{2}} , displaystyle{sum_{i , = , 1}^2} ,||{cal A}_i (X) ,-, U_i^k||_F^2 end{array}$

where $U_1^k mathrel{mathop:}= -left( Z^k + (1/rho)Y_1^k right)$ and $U_2^k mathrel{mathop:}= G - (1/rho)Y_2^k$ . From first order optimality conditions we have

$-,X^{-1} ,+, rho, {cal A}_1^dagger ( {cal A}_1 (X) - U_1^k ) ,+, rho, {cal A}_2^dagger ( {cal A}_2 (X) - U_2^k ) ;=; 0.$

Since ${cal A}_1^dagger {cal A}_1$ and ${cal A}_2^dagger {cal A}_2$ 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 X_i

and adding a quadratic penalty on the difference between

and

, $X_{i+1}$ is obtained as the minimizer of

$-{rm log,det} X ; + ; rho, displaystyle{sum_{j , = , 1}^2} left<{cal A}_j^daggerleft({cal A}_j (X_i) - U_j^k right), Xright> ; +; displaystyle{frac{mu}{2}} , ||X ,-, X_i||_F^2.$

To ensure convergence of the proximal gradient method, the parameter

has to satisfy $mu ge rho, lambda_{max}( {cal A}_1^dagger {cal A}_1 + {cal A}_2^dagger {cal A}_2)$ , where we use power iteration to compute the largest eigenvalue of the operator ${cal A}_1^dagger {cal A}_1 + {cal A}_2^dagger {cal A}_2$ .

By taking the variation with respect to

we arrive at the first order optimality condition

$mu, X ,-, X^{-1} ;= ; mu, X_i ,-, rho, displaystyle{sum_{j , = , 1}^2} {cal A}_j^daggerleft({cal A}_j (X_i) - U_j^k right) hspace{1.5cm} {rm (*)}$

$displaystyle{g_j ;=; displaystyle{frac{lambda_j}{2mu}} ,+, sqrt{left(displaystyle{frac{lambda_j}{2mu}}right)^2 ,+, displaystyle{frac{1}{mu}}}}.$

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 $X_0 mathrel{mathop:}= X^k$ , we obtain $X^{k+1}$ by repeating inner iterations until the desired accuracy is reached.

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

Purpose

Problem formulation

Alternating Minimization Algorithm (AMA)

Alternating Direction Method of Multipliers (ADMM)

Acknowledgements