The most recent methods in quantum chemical geometry optimization use the computed energy and its first derivatives with an approximate second derivative matrix. The performance of the optimization process depends highly on the choice of the coordinate system. In most cases the optimization is carried out in a complete internal coordinate system using the derivatives computed with respect to Cartesian coordinates. The computational bottlenecks for this process are the transformation of the derivatives into the internal coordinate system, the transformation of the resulting step back to Cartesian coordinates, and the evaluation of the Newton-Raphson or rational function optimization (RFO) step. The corresponding systems of linear equations occur as sequences of the form y(i)=M(i)x(i), where M-i can be regarded as a perturbation of the previous symmetric matrix Mi-1. They are normally solved via diagonalization of symmetric real matrices requiring O(N-3) operations. The current study is focused on a special approach to solving these sequential systems of linear equations using a method based on the update of the inverse of the symmetric matrix Mi. For convergence, this algorithm requires a number of O(N-2) operations with an O(N-3) factor for only the first calculation. The method is generalized to include redundant (singular) systems. The application of the algorithm to coordinate transformations in large molecular geometry optimization is discussed.