CVXOPT and portfolio optimization: puzzling issue

Question

I am trying to solve the following simple optimization problem. One seeks to find the global minimum variance portfolio, being the portfolio that minimizes variance with only one constraint : weights must sum to one. Optimization program

This problem has a well-known closed-form solution: Solution

I'm trying to reproduce the results using CVXopt in Python, and I encounter a puzzling issue.

I have a 10x10 sample covariance matrix.

• If I solve the problem with the entire 10x10 matrix, the output is incorrect : the "sum-to-one" constraint is not respected, and weights are different from the closed-form solution
• If I solve the problem with a 7x7 subsample of the same matrix, the output is correct : the "sum-to-one" constraint is respected, and weights are equivalent to the closed-form solution

Actually, with any subsample of size equal or lower than 7x7, it works. But for any subsample of size higher or equal to 8x8, it does not work anymore. I can't grasp where the problem is coming from.

Thank you for your help!

Here is a sample code

%reset -sf #Clear Environment
import numpy as np
import cvxopt
from cvxopt import matrix as dmatrix
from cvxopt.solvers import qp, options
from numpy.linalg import inv
from numpy import ones
from numpy import transpose as t
# Sample 10x10 covariance matrix 
sigmafull = np.array([[0.01449082, 0.00846992, 0.00846171, 0.00773097, 0.00878925,
        0.00748843, 0.00672341, 0.00665912, 0.0068593 , 0.00827341],
       [0.00846992, 0.00952205, 0.00766057, 0.00726647, 0.00781524,
        0.00672368, 0.00642426, 0.00609368, 0.00617965, 0.00704281],
       [0.00846171, 0.00766057, 0.00842194, 0.00700168, 0.00772423,
        0.0061137 , 0.00612574, 0.00601041, 0.00621007, 0.00712152],
       [0.00773097, 0.00726647, 0.00700168, 0.00687784, 0.00726901,
        0.00573606, 0.00567145, 0.00556391, 0.00575279, 0.00660916],
       [0.00878925, 0.00781524, 0.00772423, 0.00726901, 0.00860462,
        0.00612804, 0.0061301 , 0.00603605, 0.00630947, 0.0075281 ],
       [0.00748843, 0.00672368, 0.0061137 , 0.00573606, 0.00612804,
        0.00634431, 0.0054793 , 0.00513665, 0.00511852, 0.00575049],
       [0.00672341, 0.00642426, 0.00612574, 0.00567145, 0.0061301 ,
        0.0054793 , 0.0055722 , 0.0050824 , 0.00512499, 0.00576934],
       [0.00665912, 0.00609368, 0.00601041, 0.00556391, 0.00603605,
        0.00513665, 0.0050824 , 0.00521583, 0.00510142, 0.00576414],
       [0.0068593 , 0.00617965, 0.00621007, 0.00575279, 0.00630947,
        0.00511852, 0.00512499, 0.00510142, 0.00547566, 0.00603528],
       [0.00827341, 0.00704281, 0.00712152, 0.00660916, 0.0075281 ,
        0.00575049, 0.00576934, 0.00576414, 0.00603528, 0.00756009]])
# sigma = sigmafull[0:8,0:8] #With this subsample, output is incorrect. n=8 (and for all n>8)
sigma = sigmafull[0:7,0:7] #With this subsample, output is correct. n=7 (and for all n<7)
n=len(sigma)
sigma = dmatrix(sigma) #Formatting sigma to be a dense matrix for cvxopt
mu = dmatrix(np.zeros(n)) #We just want to minimize variance, hence vector of zeroes
#Format of the equality constraint : Ax = b
#We want the sum of x to be equal to 1
Amatrix = dmatrix(ones(n)).T #Vector of ones 
bmatrix = dmatrix(1.0) #Scalar = 1
sol = qp(sigma, mu, None, None, A=Amatrix, b=bmatrix) #No inequality constraint 
w_gmv = (inv(sigma)@ones(n))/(t(ones(n))@inv(sigma)@ones(n)) #Analytical solution which indeed does sum to 1    
print(t(np.array(sol['x'])) - w_gmv) #If Vector of zeroes -> Weights are equivalent -> OK
print(sum(np.array(sol['x']))) #If equal to 1 -> OK