I doubt there is a specialized algorithm for this. This is just an equality-constrained least-squares problem, readily solved using the Lagrange multiplier approach. I don't even think you can exploit the structure of $Q$.
The Lagrangian is
$$L(x,\lambda) = \|Kx-y\|^2+\mu\|x\|^2-\lambda^T(Qx-v)$$
The optimality conditions are
$$2K^T(Kx-y)+2\mu x - Q^T\lambda=0 \quad Qx=v$$
Solving the first equation for $x$ we get
$$x=\tfrac{1}{2}(K^TK+\mu I)^{-1} Q^T \lambda$$
Substituting this into the second equation we get
$$Q(K^TK+\mu I)^{-1}Q^T \lambda = 2v$$
Assuming that Q has full row rank, you simply solve this equation for $\lambda$ and substitute the result into the previous equation for $x$. A typical algorithm would proceed as follows:
- Compute the Cholesky factorization $R_1^TR_1=K^TK+\mu I$. For a more expensive but more accurate approach, compute the QR factorization $Q_1R_1=\begin{bmatrix} K \\ \mu^{1/2} I \end{bmatrix}$.
- Compute $\tilde{Q}=QR_1^{-1}$, and then compute the Cholesky factorization $R_2^TR_2=\tilde{Q}\tilde{Q}^T$. Similarly, you can compute the QR factorization $Q_2R_2=\tilde{Q}^T$.
- Compute $\lambda = 2R_2^{-1}R_2^{-T}v$.
- Compute $x=\tfrac{1}{2}R_1^{-1}R_1^{-T}Q^T\lambda=\tfrac{1}{2}R_1^{-1}\tilde{Q}^T\lambda$.
If you're going to exploit the structure of $Q$, it will have to be in step 2, but I'm not entirely sure how you will accomplish that. Perhaps you fully invert $R_1$ and exploit the sparsity of $Q$, but I don't see how the Toeplitz structure in particular will be helpful.
A projected conjugate gradient method might allow you to exploit the full structure of $Q$, but I'm still not certain that would be a significant driver of performance.