# Tikhonov Regularized Least Squares with Unit Simplex Constraint

First, I got this $\ell_2$ norm problem, which is easy to solve

$$w_0 = \arg\min_𝑤⁡ \| 𝑋 𝑤−𝑥_𝑐 \|_2^2+ \lambda \|𝑤\|_2^2 \tag{1}$$

where $X$ is a matrix and $\lambda$ is a fixed scalar. Then, I added the following constraints

$$𝑤∈[0,1], \qquad \thinspace \sum_𝑖 𝑤_i =1 \tag{2}$$

At first I found this optimization problem similar to fmincon function of matlab (reference here), but I keep getting poor results (significant higher error than results in Eq(1)).

I already try with other matlab build-in functions such as lsqlin, lsqnonneg but all return bad results (err = $\|X*w_o - 𝑥_𝑐 \|_2^2$) such as

Eq2-fmincon: 75596.826, Eq2-lsqlin: 386.2777, Eq2-lsnonneg:311.5511, Eq1: 0.0019292

Is there any suggestion on how to solve this problem? Thank you for your time.

My matlab code for fmincon are follow

gamma = 1;
fun = @(w) ((norm(X * w  - xc, 2))^2 + gamma * (norm(w, 2))^2);
w0  = rand(size(X, 2), 1); %w0 = w0./sum(w0);
A   = [];                   b   = [];
Aeq = ones(1, size(X, 2));  beq = 1;
lb  = zeros(size(X, 2), 1); ub  = ones(size(X, 2), 1);
w1  = fmincon(fun,w0,A,b,Aeq,beq,lb,ub) ;
w2  = lsqlin(X, xc, A, b, Aeq, beq, lb, ub);
w3  = lsqnonneg(X, xc);
w4  = pinv(X' * X + gamma * eye(size(X,2))) *X' * xc;

tmp1 = norm(X * w1 - xc, 2);     tmp2 = norm(X * w2 - xc, 2);
tmp3 = norm(X * w3 - xc, 2);     tmp4 = norm(X * w4 - xc, 2);

display(['Eq2-fmincon: ' num2str(tmp1) ', Eq2-lsqlin: ' num2str(tmp2) ', Eq2-lsnonneg:' ...
num2str(tmp3) ', Eq1: ' num2str(tmp4)]);

• Fmincon is a valid approach. What algorithm did you use? Did it converge? Why do you think there are $w$ that satisfy your constraints but give lower errors? – LinAlg Oct 24 '16 at 15:39
• Thank you. I did use fmincon but got "Solver stopped prematurely" with large error. I added matlab script to my question. – Atena Nguyen Oct 24 '16 at 15:40
• Your nonlcon-code runs fine on my random data set. It might help if you provide gradient information, and if you ensure that the starting point is not too close to $0$ or $1$. – LinAlg Oct 24 '16 at 15:49
• Thank you very much. I initiate w0 via random number so it should not be close to 0 or 1. The gradient for fmincon is not stable. You can find it here dropbox.com/s/x95cio4150jj4a0/fminconGrad.png?dl=0 – Atena Nguyen Oct 24 '16 at 15:55
• Sorry. with "provide" I meant to provide it to fmincon; see "include gradient" at mathworks.com/help/optim/ug/fmincon.html It is not clear to me what your graph shows. The gradient for your problem is a vector. – LinAlg Oct 24 '16 at 15:59

## 2 Answers

This is a regularized least-squares (RLS) problem subject to the standard $(n-1)$-simplex. We have the following quadratic program (QP)

$$\begin{array}{ll} \text{minimize} & \| \mathrm A \mathrm x - \mathrm b \|_2^2 + \lambda \| \mathrm x \|_2^2\\ \text{subject to} & 1_n^{\top} \mathrm x = 1\\ & \mathrm x \geq \mathrm 0_n\end{array}$$

which can be rewritten as follows

$$\begin{array}{ll} \text{minimize} & \mathrm x^{\top} \left( \mathrm A^{\top} \mathrm A + \lambda \mathrm I_n \right) \mathrm x - 2 \mathrm b^{\top} \mathrm A \mathrm x + \mathrm b^{\top} \mathrm b\\ \text{subject to} & 1_n^{\top} \mathrm x = 1\\ & \mathrm x \geq \mathrm 0_n\end{array}$$

Do note that if $\lambda < 0$, then the QP may be non-convex.

In MATLAB, one can use function quadprog to solve this QP.

• A thousand thanks :D, i will check the matlab code soon. – Atena Nguyen Oct 25 '16 at 11:02

You can also solve this by Projected Gradient Descent since the projection onto the Unit Simplex is known.

The Code:

vX = pinv(mA) * vB;

mAA = mA.' * mA;
mAb = mA.' * vB;

for ii = 1:numIterations
stepSize = stepSizeBase / sqrt(ii);
vG = (mAA * vX) - mAb + (2 * paramLambda * vX);
vX = vX - (stepSize * vG);
vX = ProjectSimplex(vX, simplexRadius, stopThr);
end

objVal = (0.5 * sum((mA * vX - vB) .^ 2)) + (paramLambda * sum(vX .^ 2));

disp([' ']);
disp(['Projected Gradient Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(objVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX.'), ' ]']);
disp([' ']);


The full code (Validated by CVX) is in my StackExchange Mathematics Q1983030 GitHub Repository.

You can even accelerate this using Nesterov / FISTA to get really fast and efficient algorithm.