Closest Toeplitz SDP search.

% This script finds a Toeplitz Hermitian PSD matrix that is closest to a
% given Hermitian matrix, as measured by the Frobenius norm. That is, for
% a given matrix P, it solves:
%    minimize   || Z - P ||_F
%    subject to Z >= 0
%
% Adapted from an example provided in the SeDuMi documentation. Notice
% the use of SDP mode to simplify the semidefinite constraint.

% The data. P is Hermitian, but is neither Toeplitz nor PSD.
P = [ 4,     1+2*j,     3-j       ; ...
      1-2*j, 3.5,       0.8+2.3*j ; ...
      3+j,   0.8-2.3*j, 4         ];

% Construct and solve the model
n = size( P, 1 );
cvx_begin sdp
    variable Z(n,n) hermitian toeplitz
    dual variable Q
    minimize( norm( Z - P, 'fro' ) )
    Z >= 0 : Q;
cvx_end

% Display resuls
disp( 'The original matrix, P: ' );
disp( P )
disp( 'The optimal point, Z:' );
disp( Z )
disp( 'The optimal dual variable, Q:' );
disp( Q )
disp( 'min( eig( Z ) ), min( eig( Q ) ) (both should be nonnegative, or close):' );
disp( sprintf( '   %g   %g\n', min( eig( Z ) ), min( eig( Q ) ) ) );
disp( 'The optimal value, || Z - P ||_F:' );
disp( norm( Z - P, 'fro' ) );
disp( 'Complementary slackness: Z * Q, should be near zero:' );
disp( Z * Q )
 
Calling sedumi: 20 variables, 6 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 6, order n = 6, dim = 30, blocks = 3
nnz(A) = 20 + 0, nnz(ADA) = 36, nnz(L) = 21
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            2.42E+01 0.000
  1 :  -3.11E+00 4.59E+00 0.000 0.1896 0.9000 0.9000   1.37  1  1  2.2E+00
  2 :  -1.29E+00 4.06E-01 0.000 0.0886 0.9900 0.9900   1.26  1  1  2.7E-01
  3 :  -1.45E+00 1.03E-02 0.000 0.0253 0.9900 0.9900   0.95  1  1  6.6E-03
  4 :  -1.45E+00 1.89E-03 0.000 0.1836 0.9000 0.9000   1.00  1  1  1.2E-03
  5 :  -1.45E+00 4.31E-04 0.000 0.2288 0.9000 0.9000   1.00  1  1  2.7E-04
  6 :  -1.45E+00 1.66E-05 0.000 0.0384 0.9907 0.9900   1.00  1  1  1.0E-05
  7 :  -1.45E+00 3.90E-06 0.199 0.2353 0.9000 0.9131   1.00  1  1  2.4E-06
  8 :  -1.45E+00 9.24E-07 0.000 0.2370 0.9000 0.9174   1.00  1  1  5.4E-07
  9 :  -1.45E+00 1.38E-07 0.000 0.1493 0.9062 0.9000   1.00  1  1  7.9E-08
 10 :  -1.45E+00 1.37E-08 0.475 0.0995 0.9900 0.9900   1.00  2  2  7.8E-09

iter seconds digits       c*x               b*y
 10      0.1   8.1 -1.4508035035e+00 -1.4508035159e+00
|Ax-b| =   7.8e-10, [Ay-c]_+ =   2.1E-09, |x|=  1.7e+00, |y|=  1.7e+00

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    1.000E-01    2.000E-02    
Max-norms: ||b||=1, ||c|| = 6,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 885.477.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.4508
The original matrix, P: 
   4.0000             1.0000 + 2.0000i   3.0000 - 1.0000i
   1.0000 - 2.0000i   3.5000             0.8000 + 2.3000i
   3.0000 + 1.0000i   0.8000 - 2.3000i   4.0000          

The optimal point, Z:
   4.2827             0.8079 + 1.7342i   2.5575 - 0.7938i
   0.8079 - 1.7342i   4.2827             0.8079 + 1.7342i
   2.5575 + 0.7938i   0.8079 - 1.7342i   4.2827          

The optimal dual variable, Q:
   0.3366            -0.0635 - 0.2866i  -0.3051 + 0.1422i
  -0.0635 + 0.2866i   0.2561            -0.0635 - 0.2866i
  -0.3051 - 0.1422i  -0.0635 + 0.2866i   0.3366          

min( eig( Z ) ), min( eig( Q ) ) (both should be nonnegative, or close):
   1.136e-09   4.18457e-10

The optimal value, || Z - P ||_F:
    1.4508

Complementary slackness: Z * Q, should be near zero:
   1.0e-05 *

   0.2927 - 0.1116i  -0.1502 - 0.2279i  -0.2178 + 0.2247i
   0.1450 - 0.6546i  -0.5846 - 0.0000i   0.1450 + 0.6546i
  -0.2178 - 0.2247i  -0.1502 + 0.2279i   0.2927 + 0.1116i