The problem seems quite ill posed so I will assume some level of over-determinism is there.
First, note that your vector, $\mathbf{z}$ lives on a hypersphere. Thus, it is indeed on a $(k-1)$ dimensional manifold, $\mathcal{S}^{k-1}$. Then, because you don't know the scales, we could always write the problem as linear combination of $k$ unit vectors, $\{\mathbf{v}_i\in \mathbb{R}^{k}\}$, absorbing the unknowns into the linear coefficients:
$$
\mathbf{z}= \lambda_1\mathbf{v}_1 + \lambda_2\mathbf{v}_2 + \cdots + \lambda_{k}\mathbf{v}_{k}
$$
As you said you could overspecify $k$, I will, for the moment, assume that $\mathbf{v}_{k}$ is available.
Let us write those vectors in a matrix:
$$
\mathbf{V} = \begin{bmatrix}
\mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_{k-1} & \mathbf{v}_{k}
\end{bmatrix}
$$
Since both $\mathbf{z}$ and $\{\mathbf{v}_i\}$ are unit length, $\boldsymbol{\lambda}=\{\lambda_i\}$ must be unit length: $\|\boldsymbol{\lambda}\|=1$. $\mathbf{z}$ can then be written as:
$$
\mathbf{z} = \mathbf{V} \boldsymbol{\lambda}
$$
Of course $\mathbf{z}$ and $\boldsymbol{\lambda}$ are the unknowns. If we have multiple projection sets (different subsets can reconstruct $\mathbf{z}$) then $i=1\dots L$ and $L>>k$ growing $\mathbf{V}$ in width, while keeping each column unit length. In contrast, we will have a set of $\boldsymbol{\lambda}$ vectors, for each projection set: $\boldsymbol{\Lambda}=\{\boldsymbol{\lambda}_i\}$. We can then write the above relation as:
$$
\mathbf{Z}_{k\times D}=\mathbf{V}_{k \times L}\boldsymbol{\Lambda}_{L \times D}
$$
where $\mathbf{Z}=\begin{bmatrix} \mathbf{z} & \mathbf{z} & \cdots & \mathbf{z}\end{bmatrix}$ as each of the $D$ projection sets should recover the same $\mathbf{z}$.
At this stage, I propose to pick two paths:
1. Alternating Optimization
Start with initializing $\mathbf{z}$ to the expected value. To do that we assume $\lambda_i=1/L\,\,\forall i$ and compute $\mathbf{z}_i = \mathbf{V}\boldsymbol{\lambda}_0\,\,\forall i$ (At this stage you could also use your initialization). We can then iterate as:
a. $t=0$.
b. Fix $\mathbf{z}=\mathbf{z}_t$ and solve:
$$
\boldsymbol{\Lambda}^{\star}=\arg\min_{\boldsymbol{\Lambda}} \|\mathbf{Z}_t-\mathbf{V}\boldsymbol{\Lambda} \|_F
$$
c. Fix $\boldsymbol{\Lambda}_t=\boldsymbol{\Lambda}$ and solve:
$$
\mathbf{z}^{\star}=\arg\min_{\mathbf{z}} \|\mathbf{Z}_t-\mathbf{V}\boldsymbol{\Lambda} \|_F
$$
d. $\mathbf{z}_t\gets \mathbf{z}^{\star}$, $\,\,t\gets t+1$. Iterate to (b)
2. Sparse Coding
Note that under the observation that the number of projections (in high dimensions) grow with factorial, and the desired vector can be reconstructed from a minimal subset of projections (thus sparse $\boldsymbol{\lambda}$), this can potentially (and sorry if I made a mistake already) formulated as a sparse coding problem, that tries to learn $\boldsymbol{\Lambda}$ and $\mathbf{z}$ simultaneously.
From this perspective $\mathbf{Z}$ denotes a dictionary, and $\boldsymbol{\Lambda}$ a representation that we simultaneously like to learn to recover the input, $\mathbf{V}$. The result of this will hopefully be multiple $\mathbf{z}_i$ corresponding to the columns of the learned $\mathbf{Z}$. We could for instance average those to get the most probable solution.
I will take a look at the details of this and update soon.
