Is the solution of this linear system always positive definite? Let $P\in \mathbb{R}^{n\times r}$ be a submatrix (which consists of the first $r$ columns) of an arbitrary $n\times n$ orthogonal matrix ($1 < r < n$). Let $I_n$ denote the $n\times n$ identity matrix and $E_n$ denote the $n\times n$ matrix with all entries being $1$. Let $J_n := E_n -I_n$. Consider the following linear system
$$P^T(J_n\circ PSP^T)P = I_{r},$$
where ''$\circ$'' denotes the Hadamard product.
It can be shown that when $S$ is restricted to be an $r\times r$ symmetric matrix, this linear system always has a unique solution.
Question: Is the unique solution $S$ always positive definite?

FYI: Illustration of the uniqueness of the solution.
Let $\text{vec}(\cdot)$ denote the vectorization operator by stacking the columns of a matrix into a single column  and $\text{vech}(\cdot)$ denote the vectorization operator for symmetric matrices by stacking only diagonal and lower-diagonal entries column by column. Let ''$\otimes$'' denote the Kronecker product. The linear system when restricted to the symmetric matrix cone can be rewritten as
$$[H_r(P\otimes P)^T \text{Diag}(\text{vec}(J_n))(P\otimes P) G_r] \text{vech}(S) = \text{vech}(I_r),$$
where $\text{Diag}(\text{vec}(J_n))$ is an $n^2 \times n^2$ diagonal matrix with $n^2-n$ nonzero diagonal entries $1$, $P\otimes P \in \mathbb{R}^{n^2 \times r^2}$ is of  full-rank, $H_r \in \mathbb{R}^{r(r+1)/2\times r^2}$ and $G_r \in \mathbb{R}^{r^2\times r(r+1)/2}$ are full-rank matrices such that for any $r\times r$ symmetric matrix $X$, $$\text{vech}(X) = H_r\text{vec}(X), \quad \text{vec}(X) = G_r\text{vech}(X), \quad \forall \ r\times r \ \text{symmetric matrix} \ X.$$   Notice that $n^2-n\geq r^2$ for any $1 < r < n$. Since $P\otimes P$ is the first $r^2$ columns of an orthogonal matrix, $\text{Diag}(\text{vec}(J_n)) (P\otimes P)$ is of full rank and thus we obtain the uniqueness.
 A: I think the answer is yes. Write $$P^\mathrm{t}(J_n \circ PSP^\mathrm{t})P=P^\mathrm{t}(PSP^\mathrm{t}-I_n\circ PSP^\mathrm{t})P=S-P^\mathrm{t}(I_n \circ PSP^\mathrm{t})P$$
so the solution satisfies
$$(1)\quad S=I_r+P^\mathrm{t}(I_n \circ PSP^\mathrm{t})P.$$
Suppose $S$ is symmetric and let $Sv_j=\lambda_j v_j$ for $j=1,2,\ldots,r$ be an orthonormal eigenbasis of $\mathbb{R}^r$. Then
$$I_n \circ PSP^\mathrm{t}=\sum_{ij}\lambda_j (e_i e_i^\mathrm{t})\circ (Pv_j v_j^\mathrm{t}P^\mathrm{t})=\sum_{ij}\lambda_j (e_i \circ Pv_j)(e_i \circ Pv_j)^\mathrm{t}$$
$$=\sum_{i} (e_i e_i^\mathrm{t}) \sum_j \lambda_j (P_i \cdot v_j )^2.$$
Furthermore
$$\left\langle v_i, P^\mathrm{t}(I_n \circ PSP^\mathrm{t})Pv_i\right\rangle = \left\langle Pv_i, (I_n \circ PSP^\mathrm{t})Pv_i \right\rangle\geq \|v_i\|^2\min_i \sum_j \lambda_j(P_i \cdot v_j )^2$$
which establishes that the smallest eigenvalue of $P^\mathrm{t}(I_n \circ PSP^\mathrm{t})P$ is no larger than the smallest eigenvalue of $S$. Therefore, if $S$ has a negative eigenvalue, the statement (1) leads to a contradiction. I think this means that $S$ and $S-I_r$ are positive semidefinite, which would imply that $S$ is positive definite. A similar argument would apply if $I_r$ were replaced by a different positive definite matrix.
