Projection matrices with constant diagonal Given $k<n$, does there always exist an $n\times n$ projection matrix of rank $k$ with a constant diagonal, i.e., all diagonal entries equal to $k/n$? Here, when I say projection matrix, I mean specifically orthogonal projection, so $P=P^2=P^*$
It is easy to construct examples when $k=1$, or when $k=2$ and $n$ is even, and given a Hadamard matrix of order $n$, one can use it to construct examples of projection matrices of rank k and size n for all $k<n$.  However, trying to find an example with $n=3, k=2$ seems difficult enough (the way I approached it required using a computer to solve a system of 6 quadratic equations in 6 unknowns).
Is there a construction that works for general $n,k$? Is there a non-constructive proof that such a projection exists for every $n$ and $k$?
 A: As I note in my comment on your question, corollary 5 in these notes leads to a non-constructive proof that such projections must exist for every size and rank. That said, here's a nice approach for constructing constant-diagonal projections of arbitrary size and rank.
Of course, the case of $k=n$ is simple: the identity matrix is the only choice. Similarly, if $k = 0$, the $0$ matrix is the only choice. For the remaining cases, let $v_{0},v_{1},\dots,v_{n-1}$ denote the columns of the size-$n$ DFT matrix. If $k$ is even, we can take
$$
M = \sum_{j=1}^{k/2} v_j v_j^* + \sum_{j=1}^{k/2} v_{n-j}v_{n-j}^* = 2\operatorname{Re}\left[ \sum_{j=1}^{k/2} v_j v_j^*\right],
$$
where for a matrix $A$, $A^*$ denotes its conjugate-transpose and $\operatorname{Re}[A]$ denotes the (entrywise) real part of $A$. The entries of this matrix are given by
$$
M_{pq} = \frac 2n\sum_{j=1}^{k/2} \cos(2 \pi (p-q)j).
$$
If $k$ is odd, we can take
$$
M = v_0v_0^* + \sum_{j=1}^{(k-1)/2} v_j v_j^* + \sum_{j=1}^{(k-1)/2} v_{n-j}v_{n-j}^* = v_0v_0^* + 2\operatorname{Re}\left[ \sum_{j=1}^{k/2} v_j v_j^*\right].
$$
The entries of this matrix are given by
$$
M_{pq} = \frac 1n + \frac 2n\sum_{j=1}^{k/2} \cos(2 \pi (p-q)j).
$$

Bonus: a Python script
import numpy as np
import scipy.linalg as la

def proj(n,k):
    F = la.dft(5)
    U = F[:,1:(1 + k//2)]
    U = np.asmatrix(U)
    M = 2*np.real(U @ U.H)/n
    if k%2:
        M += 1/n
    return M

print(proj(5,3))

