I would use Accelerated Gradient Descent Method.
This means you'll stay only with Matrix Multiplication (And exploiting $ {A}_{i} $ being Symmetric can help with that).
But you actually need to use "Projected Sub Gradient", namely to calculate the projection of the solution onto the Domain.
The Projection onto The Constrain
Let's define a vector $ z = \left[ {{x}_{1}}^{T}, {{x}_{2}}^{T}, \ldots, {{x}_{K}}^{T} \right]^{T} $.
The constrains is equivalent of:
$$ S * z = \boldsymbol{1}_{n}, \; S = \underset{\times K}{\left [ \underbrace{{I}_{n}, {I}_{n}, \ldots, {I}_{n}} \right ]} $$
Basically the matrix $ S $ sum over the $ i $ -th element of all $ x $ vectors.
Now, given a vector $ y \in \mathbb{R}^{nk} $, its projection onto the set $ \mathcal{S} = \left\{ x \mid S * x = \boldsymbol{1}_{n} \right\} $ is given by:
$$ \arg \min_{x \in \mathcal{S}} \frac{1}{2} \left\| x - y \right\|_{2}^{2} = y - {S}^{T} \left( S {S}^{T} \right)^{-1} \left( S y - \boldsymbol{1}_{n} \right) $$
Due to the special structure of $ S $ one could see it is equivalent to:
$$ {x}_{i} - \frac{\sum_{i = 1}^{K} {x}_{i} - \boldsymbol{1}_{n}}{K}, \: i = 1, 2, \ldots, k $$
Namely spread the deviation equally on all elements.
Here is the code:
%% Solution by Projected Gradient Descent
mX = zeros([numRows, numCols]);
for ii = 1:numIterations
for jj = 1:numCols
mX(:, jj) = mX(:, jj) - (stepSize * ((2 * tA(:, :, jj) * mX(:, jj)) + (paramLambda * mC(:, jj))));
end
mX = hProjEquality(mX);
end
objVal = 0;
for ii = 1:numCols
objVal = objVal + (mX(:, ii).' * tA(:, :, ii) * mX(:, ii)) + (paramLambda * mC(:, ii).' * mX(:, ii));
end
disp([' ']);
disp(['Projected Gradient Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(objVal)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);
The full code is in my Mathematics Q2199546 GitHub Repository (Specifically in Q2199546.m).
The code is validated using CVX.
Best Answer
Yes, from the observation that moving two $x_i$ "closer together" without changing the sum—that is, if $x_i < x_j$, then replacing $x_i$ and $x_j$ with $x_i + \epsilon$ and $x_j - \epsilon$—will always reduce $\sum x_i^2$. It's easy to see this by writing $x_i + x_j = k$ and rewriting $x_i^2 + x_j^2 = x_i^2 + (k-x_i)^2 = 2 \left(x_i - \frac{k}{2}\right)^2 + \frac{k^2}{2}$, which is minimized when $x_i = k/2$ and increases with $\left|\frac{k}{2} - x_i\right|$.
This means that the (necessarily unique) arrangement of $x_i$ that gives a minimum value of $\sum x_i^2$ cannot have any $x_i$ that is less than both its corresponding $a_i$ and some other $x_j$; otherwise, there would be room to reduce $\sum x_i^2$ by raising $x_i$ and lowering $x_j$. This forces a minimum of the form that you give.
One note: it's possible for the minimizer of $\sum x_i^2$ to give $x_i = a_i$ even when $a_i > m/n$. Consider the case $m = 1$, $n = 3$, $a_1 = 0.2$, $a_2 = 0.35$, $a_3 = 1$. "Large enough" doesn't necessarily mean greater than $m/n$.