Given $[D]$ and $\vec p$, can we solve $[D]=\begin{bmatrix}p\\q\end{bmatrix}^{\top}\begin{bmatrix}c\\1...\end{bmatrix}$ for $\vec q$ and $\vec c$? Given matrix $[D]\in \Bbb R^{m \times n}$ of spectral data composed of $n$ mixtures of two pure component spectra $\vec p$ and $\vec q\in\Bbb R^{m}$, scaled by matrix $[C]\in\Bbb R^{2 \times n}$ consisting of the concentrations of the two components in each mixture. The $m$-dimension represents the spectrum x-axis index and the $n$-dimension represents the sample mixture index. We have: $$[D]_{m\times n} = \begin{bmatrix}\vec p \\\vec q \end{bmatrix}^{\top}_{m\times2}[C]_{2\times n}$$
Imagine a set of spectra $[D]$ that have been normalized by the concentration of the second component and that we know the first pure component spectrum $\vec p$. Thus we have: $$[D]_{m\times n} = \begin{bmatrix}\vec p \\\vec q \end{bmatrix}^{\top}_{m\times2}\begin{bmatrix}\,\vec c\\1\cdots 1\end{bmatrix}_{2\times n}$$How can we obtain the least-squares estimates of $\vec q$ and $\vec c$?
 A: The notation $[D]: \mathbb R \to \mathbb R^{m \times n}$ confuses me, so I suppose you've meant simply $D \in \mathbb R^{n \times m}$. I also will treat $p, q \in \mathbb R^m$ and $c \in \mathbb R^n$ as column vectors.
Let's simplify the matrix notation:
$$
D = 
\begin{bmatrix}
p^\top\\
q^\top
\end{bmatrix}^\top
\begin{bmatrix}
c^\top\\
e^\top
\end{bmatrix} =
\begin{bmatrix}
p & q
\end{bmatrix}
\begin{bmatrix}
c^\top\\
e^\top
\end{bmatrix} = 
p c^\top + q e^\top.
$$
Here $e \in \mathbb R^n$ denotes a vector of ones. It is now obvious that the right hand side is linear both in $c$ and $q$.
The least squares solution minimizes the residual's sum of squares that is
$$
E = \sum_{i=1}^m \sum_{j=1}^n (d_{ij} - p_i c_j - q_i)^2 \to \min_{c_j, q_i}.
$$
This is a quadratic problem and the solution can be obtained from the optimality conditions:
$$
0 = \frac{\partial E}{\partial c_j} = \sum_{i=1}^m 2 (d_{ij} - p_i c_j - q_i) (-p_i) = 
-\sum_{i=1}^m p_i (d_{ij} - p_i c_j - q_i) = -2p^\top (D - pc^\top - qe^\top)\\
0 = \frac{\partial E}{\partial q_i} = -\sum_{j=1}^n 2 (d_{ij} - p_i c_j - q_i) =
-2 (D - pc^\top - qe^\top) e = 0.
$$
Transposing the first equation we have
$$
(D^\top - cp^\top - eq^\top) p = 0.
$$
Now separating unknowns from the knowns we obtain a system of linear equaions:
$$
(cp^\top + e q^\top) p = D^\top p\\
(pc^\top + qe^\top) e = D e
$$
Using property $a^\top b = b^\top a$ when the product is scalar and noting that $a^\top b$ as a scalar can be safely moved across the product we get
$$
(p^\top p) c + (e p^\top) q = D^\top p\\
(pe^\top) c + (e^\top e) q = De
$$
In matrix form this system can be written as
$$
\begin{bmatrix}
p^\top p I & ep^\top\\
pe^\top & n I
\end{bmatrix}
\begin{bmatrix}
c\\q
\end{bmatrix} =
\begin{bmatrix}
D^\top p\\De
\end{bmatrix}
$$
This system's matrix is singular, but consistent. It can be checked by multiplying with $[e^\top\; -p^\top]$ on the left. The solution is also not unique. If $c_0, q_0$ is a solution then
$$
c = c_0 + \alpha e, \quad q = q_0 - \alpha p.
$$
also would be. $c_0, q_0$ may be obtained using the pseudoinverse matrix.
Note that for arbitrary $\alpha$ the difference between $D$ and $pc^\top + qe^\top$ will be the same due to
$$
pc^\top + qe^\top = pc_0^\top + \alpha pe^\top + q_0 e^\top - \alpha pe^\top = 
pc_0^\top + q_0 e^\top.
$$
