J. M. es correcto que la dificultad es que C no necesariamente tienen columna completa en la clasificación, la causa de su aumentada del sistema a ser singular. En el método estándar, en este caso sería para regularizar su problema, es decir, reemplazar la restricción $C^T q + c_0 = 0$ $C^T q + \delta r + c_0 = 0$ para un valor pequeño de $\delta > 0$ (por ejemplo, entre el$10^{-6}$$10^{-8}$) y una nueva variable $r$ que representan el valor residual de sus limitaciones. La restricción de la matriz es ahora
$$
\begin{bmatrix}
C^T & \delta I
\end{bmatrix}
$$
que siempre tiene una fila completa de rango. Al mismo tiempo, usted querrá agregar un plazo $\delta^2 \|r\|^2$ a su objetivo con el fin de minimizar este residual.
Si todo es nonsingular, usted tiene varias opciones para resolver el problema. Uno es simétrica indefinida de la factorización de la aumentada del sistema (el que ustedes llaman el "normal ecuaciones"). Esto se puede hacer de manera muy eficiente y en paralelo, si el sistema es tan grande. En Matlab, mirar en la lbl función. En Fortran, mirar en la HSL biblioteca: http://www.hsl.rl.ac.uk/academic.html (los paquetes relevantes son MA27, MA47, MA57 y MA86).
Con respecto a métodos iterativos, aparte de MINRES, usted puede también intentar SYMMLQ si el sistema es nonsingular. Ver a Mike Saunders' sitio web de enlaces a esos códigos: http://www.stanford.edu/group/SOL/software.html
Hay una alternativa si usted tiene un medio para identificar de manera eficiente una solución particular $\bar{q}_0$ $C^T q + c_0 = 0$ (por factorización o el uso de alguna estructura especial del problema). Usted puede "cambiar" su sistema lineal mediante la definición de $\bar{q} := q - \bar{q}_0$:
$$
\begin{bmatrix}
M & C \\
C^T & 0
\end{bmatrix}
\begin{bmatrix}
\bar{q} \\ \lambda
\end{bmatrix}
=
\begin{bmatrix}
M(q_0 - \bar{q}_0) \\ 0
\end{bmatrix}.
$$
Estas son las condiciones de optimalidad de la (sin restricciones) lineal de mínimos cuadrados problema
$$
\min_{\lambda} \ \tfrac{1}{2} \| C \lambda - M(q_0 - \bar{q}_0) \|_{M^{-1}}^2.
$$
AHORA usted puede utilizar el conjugado de gradientes, o mejor aún, LSQR (todavía de Michael Saunders' de la página). LSQR es una variante de la CG que será más precisa si $C$ es mal condicionado. Una vez que haya encontrado $\lambda$, usted puede recuperar los $\bar{q} = q_0 - \bar{q}_0 - M^{-1} C \lambda$ y su solución es $q = \bar{q} + \bar{q}_0$.
Para obtener más general $M$ (no en diagonal), hay una buena alternativa híbrida. Consiste en la interpretación de las desplazado sistemas lineales de las condiciones de optimalidad de la cuadrática programa
$$
\min_{\bar{q}} \ -(q_0 - \bar{q}_0)^T M \bar{q} + \tfrac{1}{2} \bar{p}^T M \bar{q} \quad
\text{s.t.} \ C^T \bar{q} = 0.
$$
El método es un híbrido entre una factorización y CG. Es necesario factorizar la matriz
$$
\begin{bmatrix}
I & C \\
C^T & 0
\end{bmatrix}
$$
y se aplica CG a $M$. En cada iteración, la dirección de búsqueda se proyecta en el nullspace de $C^T$, asegurando que todos los recorre en satisfacer $C^T \bar{q} = 0$. La principal referencia es http://epubs.siam.org/sisc/resource/1/sjoce3/v23/i4/p1376_s1
Claramente, para la diagonal $M$, no hay ninguna ventaja en comparación con un directo de la factorización de su sistema original. Pero en muchas aplicaciones, la sustitución de $M$ $I$ da lugar a una matriz dispersa. También se $M$ no necesita ser positiva definida, sino que sólo positiva definida en el nullspace de $C^T$.
