% 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, 14 equality constraints
------------------------------------------------------------
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 = 14, order n = 6, dim = 30, blocks = 3
nnz(A) = 28 + 0, nnz(ADA) = 196, nnz(L) = 105
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            6.85E+00 0.000
  1 :   2.77E+00 1.30E+00 0.000 0.1896 0.9000 0.9000   0.92  1  1  9.8E-01
  2 :   1.44E+00 2.38E-01 0.000 0.1834 0.9000 0.9000   1.22  1  1  2.6E-01
  3 :   1.45E+00 7.18E-03 0.000 0.0301 0.9900 0.9900   1.02  1  1  8.2E-03
  4 :   1.45E+00 1.63E-04 0.000 0.0227 0.9900 0.9900   1.00  1  1  1.9E-04
  5 :   1.45E+00 3.25E-05 0.000 0.1988 0.9000 0.9000   1.00  1  1  3.7E-05
  6 :   1.45E+00 1.20E-06 0.000 0.0371 0.9900 0.8832   1.00  1  1  1.2E-06
  7 :   1.45E+00 1.54E-07 0.201 0.1279 0.9146 0.9000   1.00  1  1  2.5E-07
  8 :   1.45E+00 1.39E-08 0.000 0.0900 0.9090 0.9000   1.00  1  1  4.3E-08
  9 :   1.45E+00 1.13E-10 0.304 0.0082 0.9990 0.9990   1.00  1  1  4.2E-10

iter seconds digits       c*x               b*y
  9      0.1   Inf  1.4508035177e+00  1.4508035180e+00
|Ax-b| =   5.8e-11, [Ay-c]_+ =   3.8E-10, |x|=  9.4e+00, |y|=  1.4e+00

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    6.000E-02    0.000E+00    
Max-norms: ||b||=4, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.065.
------------------------------------------------------------
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.5574 - 0.7938i
   0.8079 - 1.7342i   4.2827             0.8079 + 1.7342i
   2.5574 + 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.85125e-09   -3.81885e-10

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

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

   0.3750 - 0.0708i  -0.1311 - 0.3062i  -0.3103 + 0.2227i
   0.1860 - 0.8396i  -0.7503 + 0.0000i   0.1860 + 0.8396i
  -0.3103 - 0.2227i  -0.1311 + 0.3062i   0.3750 + 0.0708i