# An Efficient Algorithm for the Separable Nonlinear Least Squares Problem
## 1 Introduction  
Many applications [1; 2; 3] lead to the need to find a solution \((y^{*},z^{*})\) of the separable nonlinear overdetermined least squares problem  
\[min_{y,z}\|F(y,z)\|\equiv min_{y,z}\|A(y)z+b(y)\|, \tag{1}\]  
where the full-rank matrix \(A(y)\in\mathbf{R}^{(N+\ell)\times N}\) and the vector \(b(y)\in\mathbf{R}^{N+\ell}\) are twice Lipschitz continuously differentiable functions of \(y\in\mathbf{R}^{n}\), \(z\in\mathbf{R}^{N}\) and \(\ell\geq n\). Here, we are using the Euclidean norm. The classic Golub-Pereyra variable projection method [4] replaces the problem (1) by the simpler least squares problem  
\[min_{y}\|f(y)\|\equiv min_{y}\|[I-A(y)A^{+}(y)]b(y)\|. \tag{2}\]  
Once \(y^{*}\) has been obtained by solving (2), \(z^{*}\) can be found by solving the resulting linear least squares problem  
\[min_{z}\|A(y^{*})z+b(y^{*})\|; \tag{3}\]  
formally,  
\[z^{*}=-A^{+}(y^{*})b(y^{*}), \tag{4}\]  
where \(A^{+}(y^{*})\) is the pseudo-inverse of \(A(y^{*})\). An alternative method that reduces (1) to a smaller least squares problem \(min_{y}\|f(y)\|\) was proposed in [5; 6] for the case \(\ell=n\), and the associated iterativetechnique was then applied to the case \(\ell>n\) in [7] but without a complete analytical justification. The reduced system is  
\[min_{y}\|f(y)\|\equiv min_{y}\|C^{T}(y)b(y)\|, \tag{5}\]  
where \(C(y)\) is an \((N+\ell)\times N\) matrix whose columns form an orthonormal basis for \(Null(A^{T}(y))\). In [8], we presented a quadratically convergent Gauss-Newton type iterative method, incorporating second derivative terms, to solve this problem, and provided a complete theoretical justification for our technique. In particular, we showed in [8] that our Gauss-Newton type method, and also the standard Gauss-Newton method [9; 10; 11], which omits the second derivative terms and is not quadratically convergent, are both invariant with respect to the specific basis matrix \(C(y)\) that is used at any particular point in the iteration. This makes it possible to substitute, at each point of the iteration, any easily computed local orthonormal basis for the nullspace of \(A^{T}(y)\). This observation also makes it possible to compute both the first and the second derivatives involved in the computation very cheaply. Many instances of (1), such as those arising from the discretization of linear differential equations with only a few nonlinear parameters, involve \(N\rightarrow\infty\), while \(\ell\) and \(n\) are fixed and small. In this case, the main computational cost in each step of our quadratically convergent Gauss-Newton type method of [8] is the QR factorization of \(A(y)\), costing approximately \(4N^{3}/3\) flops [12] per iteration. In this note, we show how the relevant computations can instead be performed by using one LU factorization of \(A(y)\) per iteration, so that the computational cost of the computation is essentially halved to \(2N^{3}/3\) flops [12]. We note that, in the case of discretizations of differential equations, the matrix \(A(y)\) is typically very sparse and strongly patterned. Since LU factorization more easily exploits and preserves such sparsity and patterns than does QR factorization, this technique also opens up opportunities for further reductions in the computational cost in such applications. In a previous paper [3], we presented the use of LU factorizations in the context of solving separable systems of nonlinear equations and zero-residual separable nonlinear least squares problems. We mentioned there that extending that approach to nonzero residual problems remained an open question. The current paper resolves that question, while also improving on the technique of [3] in other ways. In particular, our new method does not involve the singular value decomposition, and in contrast to the use of the standard Gauss-Newton method in [3], it retains quadratic convergence even in the case of nonzero residual problems. The latter is achieved by retaining the second derivative terms, and we show here how those terms can be computed very efficiently.