Norm inequality for positive semidefinite matrices Suppose $A, B\ge0$ are positive semidefinite matrices on the complex field, is it true that
$$\Vert A^2 +B^2 \Vert \le \Vert A + B \Vert^2,$$ for the spectral norm? I have tried numeric tests and the inequality seems valid. But I could not prove it. Furthermore I want to know if this inequality can be extended to finite number of PSD matrices: $$\Vert \sum_i B_i^2 \Vert\le \Vert \sum_i B_i \Vert^2,$$ for $B_i \ge 0$?
 A: Since $-(A+B)\le A-B\le A+B$, we have $\|A-B\|_2\le\|A+B\|_2$. It follows that
\begin{aligned}
\|A^2+B^2\|_2
&=\frac12\left\|(A-B)^2+(A+B)^2\right\|_2\\
&\le\frac12\left(\left\|(A-B)^2\right\|_2+\left\|(A+B)^2\right\|_2\right)\\
&=\frac12\left(\|A-B\|_2^2+\|A+B\|_2^2\right)\\
&\le\frac12\left(\|A+B\|_2^2+\|A+B\|_2^2\right)\\
&=\|A+B\|_2^2.
\end{aligned}
By Euler's four-square identity and its generalisation to sums of eight squares, I think the above proof can be generalised to the case when the LHS of the inequality is a sum of eight squares, but I am not sure if it can be generalised further.
A: This sketch of proof for the generic case is inspired by @user1551's answer.
The key is to express the sum of squares as an average of complete squares. For summations involving $n$ matrices, we wish to find $s_{ji}=\pm1$, s.t.:
$$\sum_i^n B_i^2 =\frac{1}{m} \sum_{j=1}^m (\sum_{i=1}^n s_{ji} B_i)^2.$$
Then
$$
\begin{aligned}
\Vert \sum_i^n B_i^2 \Vert &=\frac{1}{m} \Vert \sum_{j=1}^m (\sum_{i=1}^n s_{ji} B_i)^2\Vert \le \frac{1}{m} \sum_{j=1}^m  \Vert (\sum_i s_{ji} B_i)^2\Vert  \\
&=\frac{1}{m} \sum_{j=1}^m  \Vert \sum_i s_{ji} B_i\Vert^2 \le \frac{1}{m} \sum_{j=1}^m  \Vert \sum_i B_i\Vert^2 =  \Vert \sum_i B_i\Vert^2.
\end{aligned}
$$
The proof equates to find the $m\times n$ matrix $S$ composed of $\pm 1$s and whose column vectors are mutually orthogonal. Since $m$ can be arbitrarily large, this basis set is always possible. Here I give an explicit construction for $m=2^{n-1}$. For
$ S =( \vec s_1, \vec s_2, \cdots \vec s_n) $, we choose $\vec s_1=(1,\ldots,1)^T$, $\vec s_2$ to flip the sign of the second half of $\vec s_1$, $\vec s_2$ to flip the sign of second and fourth quarter of $\vec s_1$ and so on till eventually $\vec s_n =(1,-1,\cdots 1,-1)^T$.
