bug report: `lsq_classopath` function yields inferior solution

[Hua-Zhou]

Reported by Peng Zeng @ Auburn University

To reproduce the problem:

x = [ones(64, 1), (fullfact([2 2 2 2 2 2]) - 1.5) / 4.0];
b = [0 20 -18 -2 2 5 9]';
y = x * b;
 
Cmat = [0 1 1 1 0 0 0; 0 0 1 0 1 0 0];
dvec = [0; -20];
Emat = [0 0 0 1 1 0 0];
fvec = 0;
penidx = [false; true(6, 1)];   % intercept unpenalized
 
[ldapath, bpath] = lsq_classopath(x, y, -Cmat, -dvec, Emat, fvec, ...
                        'qp_solver', 'matlab', 'init_method', 'qp', ...
                        'penidx', penidx, 'epsilon', 1e-8)
 
i = 5;
0.5 * sum((y - x * bpath(:, i)).^2) + ldapath(i) * sum(abs(bpath(penidx, i))) 
 
% set estimate of intercept as zero
bpath2 = bpath; bpath2(1, :) = 0.0;
0.5 * sum((y - x * bpath2(:, i)).^2) + ldapath(i) * sum(abs(bpath2(penidx, i))) 

At the 5th tuning parameter value, setting intercept estimate to 0 gives a better objective value (214) than the solution from lsq_classopath (215.7578). This can also be confirmed by the estimate by lsq_constrsparsereg at the fixed tuning parameter

[beta5, stat] = lsq_constrsparsereg(x, y, ldapath(i), 'A', -Cmat, 'b', -dvec, ...
    'Aeq', Emat, 'beq', fvec, 'penidx', penidx, 'method', 'qp')
0.5 * sum((y - x * beta5).^2) + ldapath(i) * sum(abs(beta5(penidx)))