CCAMA – A customized alternating minimization algorithm for structured covariance completion problems
PurposeThis 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 formulationConsider a linear timeinvariant system where is the state vector, is the output, is a stationary zeromean 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 steadystate covariance matrix of the state vector if and only if the linear equation is solvable for . Here, is the expectation operator, represents the crosscorrelation of the state and the input , and denotes the complex conjugate transpose. The algebraic relation between secondorder statistics of the state and forcing can be used to explain partially known sampled secondorder statistics using stochasticallydriven 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 lowcomplexity 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 secondorder 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 lowcomplexity structures for stochastic forcing and facilitate construction of a particular class of lowpass filters that generate coloredintime 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 wellsuited 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 stepsize . The minimization step amounts to a matrix inversion and the minimization step amounts to application of the softthresholding 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 stepsize selection rules. Typically, either a constant stepsize is selected or the stepsize 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 righthandside 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. AcknowledgementsThis project is supported by the
