% Boyd & Vandenberghe, "Convex Optimization"
% Joëlle Skaf - 08/29/05
%
% Solves an extension of the classical Markovitz portfolio optimization
% problem:      minimize    x'Sx
%                   s.t.    p_'*x >= r_min
%                           1'*x = 1,   x >= 0
%                           sum_{i=1}^{0.1*n}x[i] <= alpha
% where p_ and S are the mean and covariance matrix of the price range
% vector p, x[i] is the ith greatest component in x.
% The last constraint can be replaced by this equivalent set of constraints
%                           r*t + sum(u) <= alpha
%                           t*1 + u >= x
%                           u >= 0

% Input data
randn('state',0);
n = 25;
p_mean = randn(n,1);
temp = randn(n);
sig = temp'*temp;
r = floor(0.1*n);
alpha = 0.8;
r_min = 1;

% original formulation
fprintf(1,'Computing the optimal Markovitz portfolio: \n');
fprintf(1,'# using the original formulation ... ');

cvx_begin
    variable x1(n)
    minimize ( quad_form(x1,sig) )
    p_mean'*x1 >= r_min;
    ones(1,n)*x1 == 1;
    x1 >= 0;
    sum_largest(x1,r) <= alpha;
cvx_end

fprintf(1,'Done! \n');
opt1 = cvx_optval;

% equivalent formulation
fprintf(1,'# using an equivalent formulation by replacing the diversification\n');
fprintf(1,'  constraint by an equivalent set of linear constraints...');

cvx_begin
    variables x2(n) u(n) t(1)
    minimize ( quad_form(x2,sig) )
    p_mean'*x2 >= r_min;
    sum(x2) == 1;
    x2 >= 0;
    r*t + sum(u) <= alpha;
    t*ones(n,1) + u >= x2;
    u >= 0;
cvx_end

fprintf(1,'Done! \n');
opt2 = cvx_optval;

% Displaying results
disp('------------------------------------------------------------------------');
disp('The optimal portfolios obtained from the original problem formulation and');
disp('from the equivalent formulation are respectively: ');
disp([x1 x2])
disp('They are equal as expected!');
Computing the optimal Markovitz portfolio: 
# using the original formulation ...  
Calling sedumi: 105 variables, 54 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
Split 1 free variables
eqs m = 54, order n = 82, dim = 107, blocks = 2
nnz(A) = 556 + 0, nnz(ADA) = 2212, nnz(L) = 1133
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            6.96E+00 0.000
  1 :  -1.70E+01 2.69E+00 0.000 0.3861 0.9000 0.9000   1.39  1  1  1.5E+01
  2 :  -4.56E+00 1.66E+00 0.000 0.6180 0.9000 0.9000   3.83  1  1  4.0E+00
  3 :  -1.20E-01 7.69E-01 0.000 0.4631 0.9000 0.9000   3.68  1  1  8.3E-01
  4 :   5.80E-01 1.61E-01 0.000 0.2098 0.9000 0.9000   1.76  1  1  1.3E-01
  5 :   6.96E-01 4.80E-02 0.000 0.2975 0.9000 0.9000   1.18  1  1  3.5E-02
  6 :   7.35E-01 1.37E-02 0.000 0.2851 0.9000 0.9000   1.06  1  1  9.8E-03
  7 :   7.47E-01 2.84E-03 0.000 0.2076 0.9000 0.9000   1.02  1  1  2.0E-03
  8 :   7.49E-01 6.06E-04 0.000 0.2131 0.9017 0.9000   1.00  1  1  4.2E-04
  9 :   7.50E-01 1.34E-04 0.000 0.2207 0.9049 0.9000   1.00  1  1  8.6E-05
 10 :   7.50E-01 3.19E-05 0.168 0.2384 0.9064 0.9000   1.00  1  1  1.9E-05
 11 :   7.50E-01 7.31E-06 0.105 0.2296 0.9041 0.9000   1.00  1  1  4.0E-06
 12 :   7.50E-01 1.50E-06 0.047 0.2044 0.9000 0.9019   1.00  1  1  8.5E-07
 13 :   7.50E-01 3.02E-07 0.093 0.2019 0.9000 0.9048   1.00  1  1  1.8E-07
 14 :   7.50E-01 6.58E-08 0.163 0.2179 0.9000 0.9056   1.00  2  2  4.3E-08
 15 :   7.50E-01 1.45E-08 0.105 0.2202 0.9000 0.9048   1.00  2  2  1.0E-08

iter seconds digits       c*x               b*y
 15      0.1   7.7  7.5016966445e-01  7.5016964760e-01
|Ax-b| =   1.1e-08, [Ay-c]_+ =   3.1E-09, |x|=  1.1e+00, |y|=  8.9e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    1.000E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 2.247311e+01,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1084.86.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.75017
Done! 
# using an equivalent formulation by replacing the diversification
  constraint by an equivalent set of linear constraints... 
Calling sedumi: 105 variables, 54 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
Split 1 free variables
eqs m = 54, order n = 82, dim = 107, blocks = 2
nnz(A) = 556 + 0, nnz(ADA) = 2212, nnz(L) = 1133
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            6.96E+00 0.000
  1 :  -1.70E+01 2.69E+00 0.000 0.3861 0.9000 0.9000   1.39  1  1  1.5E+01
  2 :  -4.56E+00 1.66E+00 0.000 0.6180 0.9000 0.9000   3.83  1  1  4.0E+00
  3 :  -1.20E-01 7.69E-01 0.000 0.4631 0.9000 0.9000   3.68  1  1  8.3E-01
  4 :   5.80E-01 1.61E-01 0.000 0.2098 0.9000 0.9000   1.76  1  1  1.3E-01
  5 :   6.96E-01 4.80E-02 0.000 0.2975 0.9000 0.9000   1.18  1  1  3.5E-02
  6 :   7.35E-01 1.37E-02 0.000 0.2851 0.9000 0.9000   1.06  1  1  9.8E-03
  7 :   7.47E-01 2.84E-03 0.000 0.2076 0.9000 0.9000   1.02  1  1  2.0E-03
  8 :   7.49E-01 6.06E-04 0.000 0.2131 0.9017 0.9000   1.00  1  1  4.2E-04
  9 :   7.50E-01 1.34E-04 0.000 0.2207 0.9049 0.9000   1.00  1  1  8.6E-05
 10 :   7.50E-01 3.19E-05 0.168 0.2384 0.9064 0.9000   1.00  1  1  1.9E-05
 11 :   7.50E-01 7.31E-06 0.105 0.2296 0.9041 0.9000   1.00  1  1  4.0E-06
 12 :   7.50E-01 1.50E-06 0.047 0.2044 0.9000 0.9019   1.00  1  1  8.5E-07
 13 :   7.50E-01 3.02E-07 0.093 0.2019 0.9000 0.9048   1.00  1  1  1.8E-07
 14 :   7.50E-01 6.58E-08 0.163 0.2178 0.9000 0.9056   1.00  2  2  4.3E-08
 15 :   7.50E-01 1.45E-08 0.105 0.2202 0.9000 0.9048   1.00  2  2  1.0E-08

iter seconds digits       c*x               b*y
 15      0.1   7.7  7.5016966443e-01  7.5016964760e-01
|Ax-b| =   1.1e-08, [Ay-c]_+ =   3.1E-09, |x|=  1.1e+00, |y|=  8.9e+00

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    1.100E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 2.247311e+01,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1084.84.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.75017
Done! 
------------------------------------------------------------------------
The optimal portfolios obtained from the original problem formulation and
from the equivalent formulation are respectively: 
    0.0000    0.0000
    0.0000    0.0000
    0.1342    0.1342
    0.0000    0.0000
    0.0000    0.0000
    0.1177    0.1177
    0.1134    0.1134
    0.0123    0.0123
    0.0904    0.0904
    0.0256    0.0256
    0.0451    0.0451
    0.0437    0.0437
    0.0000    0.0000
    0.1435    0.1435
    0.0000    0.0000
    0.0086    0.0086
    0.1177    0.1177
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0000    0.0000
    0.0313    0.0313
    0.1164    0.1164
    0.0000    0.0000

They are equal as expected!