机器学习中最常见的损失函数就是最小平方误差，所以我们来看看 LSE 有什么高效的解法。
LSE中，我们的优化目标常常是
$f(x) = \frac{1}{2}\sum_{j=1}^mr^2_j(x)$
令 $r(x) = (r_1(x),r_2(x),...,r_m(x))^T$，则上面的函数可以简写为 $f(x) = \|r(x)\|^2$ 。令
$J(x) = \begin{bmatrix}\nabla r_1(x)^T\\\nabla r_2(x)^T\\.\\.\\.\\\nabla r_m(x)^T\end{bmatrix}$
则
$$\begin{align}
\nabla f(x) &= \sum_{j=1}^m r_j(x)\nabla r_j(x) = J(x)^Tr(x)\\
\nabla^2 f(x) &= \sum_{j=1}^m \nabla r_j(x)\nabla r_j(x)^T + \sum_{j=1}^m r_j(x)\nabla^2 r_j(x)\\
&= J(x)^TJ(x) + \sum_{j=1}^m r_j(x)\nabla^2 r_j(x)
\end{align}$$
当我们的解接近最优解时，$r_j(x)$ 或 $\nabla^2 r_j(x)$ 很小，这时第二项可以基本忽略。
Linear Least-square
首先看一个简单的情况，我们要拟合的函数也是线性函数。这种情况下目标函数可以写成
$f(x) = \frac{1}{2} \|Jx-y\|^2$
此时一阶导和二阶导分别为：
$$\nabla f(x) = J^T(Jx-y)
\nabla^2 f(x) = J^TJ$$
求解其实只要计算 $J^TJx^* = J^Ty$ 。主要的问题在于如何计算。下面介绍3种方法计算。
Cholesky factorization
思路非常清晰：
- 计算 $J^TJ$ 和 $J^Ty$
- 对 $J^TJ$ 进行 Cholesky factorization
- 求解 $J^TX = J^Ty, Jx^* = X$ 得到 $x^*$
但是如果 $J$ 的条件数非常大，$J^TJ$ 会将其进一步放大。这种情况下 Cholesky factorization 的效果会很差，但是在 m 远大于 n 的情况下，该算法会有较好的效果
QR factorization
易知， 当 $Q$正交时， $\|Jx-y\| = \|Q(Jx-y)\|$ 根据 QR分解，有
$J\Pi = Q\begin{bmatrix}R\\0\end{bmatrix} = \begin{bmatrix}Q_1&Q_2\end{bmatrix}\begin{bmatrix}R\\0\end{bmatrix} = Q_1R$
其中， $\Pi$ 为 n n 维的 permutation 矩阵；$Q$ 为正交矩阵； $Q_1$ 是 $Q$ 的前 n 列，$Q_2$ 是 $Q$ 的剩下的 m-n 列； $R$ 为 nn 的正定上三角矩阵。所以
$$\begin{align}
\|Jx-y\|_ 2^2 &= \Bigg\|\begin{bmatrix}Q_1^T\\Q_2^T\end{bmatrix}(J\Pi\Pi^Tx-y)\Bigg\|\\
& =\Bigg\|\begin{bmatrix}R\\0\end{bmatrix}\Pi^Tx - \begin{bmatrix}Q_1^Ty\\Q_2^Tys\end{bmatrix}\Bigg\|\\
& =\|R(\Pi^Tx) - Q_1^Ty\|^2_2 + \|Q^T_2y\|^2_2
\end{align}$$
最后我们需要解 $x^* = \Pi R^{-1}Q_1^Ty$ 。实际一般解 $Rz = Q_1^Ty$ ，然后令 $x^* = \Pi z$ 。它的优点在于避开了 Cholesky factorization 中的坑，$J$ 的条件数很大时能拿到相对较好的结果。
SVD
当我们需要更为健壮的结果时，就要用到SVD了。对 $J$ 做SVD分解，可得：
$J = U\begin{bmatrix}S\\0\end{bmatrix}V^T = \begin{bmatrix}U_1&U_2\end{bmatrix}\begin{bmatrix}S\\0\end{bmatrix}V^T = U_1SV^T$
这种情况下
$$\begin{align}
\|Jx-y\| &= \Bigg\|\begin{bmatrix}S\\0\end{bmatrix}(V^Tx)-\begin{bmatrix}U_1^T\\U_2^T\end{bmatrix}y\Bigg\|^2\\
& = \|S(V^Tx)-U_1^Ty\|^2 + \|U_2^Ty\|^2
\end{align}$$
可以得到 $x^* = VS^{-1}U_1^Ty = \sum_{i=1}^n\dfrac{u_i^Ty}{\sigma_i}v_i$ 。当然，有些 $\sigma_i$ 可能是0， 这时候可以将公式改写为 $x^* = \sum_{\sigma_i\neq0}\dfrac{u_i^Ty}{\sigma_i}v_i + \sum_{\sigma_i=0}{\tau_i}v_i$ ，一般令 $\tau_i=0$ 。这是3个方法中最为健壮的。毕竟SVD分解的计算量摆在那，不好谁用啊。
当然，对于大数据量来说，还是CG什么的迭代式的更好点。
NonLinear Least-square
讲完线性的，现在来看非线性的。
Gauss-Newton method
我们可以用 $J_k^TJ_k$ 来近似二阶导，然后和牛顿法类似，解
$J_k^TJ_kp_k^{GN} = -J_k^Tr_k$
得到下降方向 $p_k^{GN}$ 。
当 $J$ 满秩， $\nabla f$ 不为0时，有
$(p_k^{GN})^T\nabla f = (p_k^{GN})^TJ_k^Tr_k = -(p_k^{GN})^TJ_kJ_k^Tp_k^{GN} = -\|J_k^Tp_k^{GN}\|^2\le0$
所以这必然是下降方向；而且很多情况下都可以认为 $B_k\approx J_k^TJ_k$ 。而且这和上面讲的线性方程的函数类似，所以上面的方法可以直接用到这里来。
当然这是 $J$ 满秩的情况。 $J$ 不满秩就比较麻烦了。这时 $p_k^{GN}$ 会有无穷多个解，上面的条件可能不再满足，因此算法的收敛性可能无法得到保证。
收敛速度方面，可以推出
$\|x_k+p_k^{GN}-x^* \| \approx \|[J^TJ(x^* )]^{-1}H(x^* )\|\|x_k-x^* \| + O(\|x_k-x^* \|^2)$
所以收敛速度方面是很快的。
在遇到大数据量的情况下，可以用 inexact newton method ，不过需要把 $B_k$ 的计算替换为 $J_k^TJ_k$
The Levenberg-marquardt method
和上面的方法类似，Levenberg-marquardt 也用了一样的方法去近似 Hessian 矩阵，但是用的搜索是 trust-region 。这样避免了 $J$ 的条件数过高带来的问题。
每次迭代，我们需要求解
$\min_p \frac{1}{2}\|J_kp+r_k\|^2\quad s.t. \|p\|\le\Delta_k$
如果 Gauss-Newton 解出的 $p^{GN}_k$ 满足 $p^{GN}_k \le \Delta_k$ ，直接用就可以了。否则根据KKT条件，需要求解
$$(J^TJ+\lambda I)p=-J^Tr\\
p = \Delta$$
可以用 trust-region 中解决子问题所用的迭代式算法去解。不过在这种场景中有更高效的 Cholesky factorization 的方法。
令 $Q_\lambda, R_\lambda$ 为正交矩阵和上三角矩阵，其中 $R_\lambda^TR_\lambda = (J^J+\lambda I)$ ，且
$$\label{eq}
\begin{bmatrix}R_\lambda\\0\end{bmatrix} = Q_\lambda \begin{bmatrix}J\\\sqrt{\lambda} I\end{bmatrix}$$
为了减小计算量，我们先对 $J$ 做QR分解：
$J = Q\begin{bmatrix}R\\0\end{bmatrix}$
对上式做一点变形，可得:
$\begin{bmatrix}R\\0\\\sqrt{\lambda}I\end{bmatrix} = \begin{bmatrix}Q^T&\\&I\end{bmatrix} \begin{bmatrix}J\\\sqrt{\lambda} I\end{bmatrix}$
对上式的左边做行变换，使其变为上三角矩阵：
$\bar{Q}_\lambda^T\begin{bmatrix}R\\0\\\sqrt{\lambda}I\end{bmatrix} = \begin{bmatrix}R_\lambda\\0\\0\end{bmatrix}$
令 $Q_\lambda = \begin{bmatrix}Q^T&\\&I\end{bmatrix}\bar{Q}_\lambda^T$
然后带到\ref{eq}就可以了。当然，对于条件数太高的情况，可以进行 scaling，这和 trust-region 基本一样，就不多说了。
大数据量和上面一样，上CG
Large-residual problem
在这种问题中，由于残差很大，所以用 $J_k^TJ_k$ 近似二阶导的效果很差，收敛速度大概为 linear，远远差于 quasi-newton 方法。但是在早期迭代时，上面方法的速度又优于 quasi-newton 法。
这种情况下我们有2种处理办法。
一是用启发式的方法确定是否使用 quasi-newton 还是 Gauss-Newton。 如果使用 Gauss-Newton 带来的损失函数减少量较大(例如减少10%)，就使用 Gauss-Newton ，并令 $B_k=J^T_kJ_k$ ，否则继续使用 $B_k$ 。但是更新 $B_k$ 时需要用类似 BFGS 的方法
二是令 $B_k = J^T_kJ_k + S_k$ 每次迭代时更新 $S_k$，更新规则为：
$S_{k+1} = S_k + \dfrac{(\hat{y}-S_ks)y^T + y(\hat{y}-S_ks)^T}{y^Ts} - \dfrac{(\hat{y}-S_ks)^Ts}{(y^Ts)^2}yy^T$
其中 $s = x_{k+1}-x_k, y = J_{k+1}^Tr_{k+1}-J_{k}^Tr_{k}, \hat{y} = J_{k+1}^Tr_{k+1}-J_{k}^Tr_{k+1}$
但是该算法不保证在 $S_k$ 会趋向于0。所以可以在 $S_k$ 之前乘上一个系数 $\tau_k = \min(1, \dfrac{|S^T\hat{y}|}{|S^TS_kS|})$
Orthogonal distance regression
在上面的问题中，我们的假设都是 t 是精确的，而 y 是不精确的。但是如果 x 也不精确，该如何处理？
这种情况下我们假设 $\epsilon_k ,\delta_k$ 分别表示 t 和 y 的误差。则实际的函数可以表示为
$y_j = \phi(x;t_j+\delta_j)+\epsilon_j$
则我们的目标为
$\min_{x,\epsilon_j,\delta_j} \frac{1}{2}\sum_{j=1}^m w_j^2\epsilon^2_j + d_j^2\delta_j^2$
上式可以重写为
$\min_{x,\delta} \frac{1}{2} \sum_{j=1}^m w_j^2[y_j-\phi(x;t_i+\delta_j)]^2 + d_j^2\delta_j^2$
令
$r_j(x,\delta) = \begin{cases}w_j^2[y_j-\phi(x;t_i+\delta_j)]^2&j=1,2...m\\d_{j-m}^2\delta_{j-m}^2&j=m+1,m+2...2m\end{cases}$
则可以化简为
$\frac{1}{2}\sum_{j=1}^{2m}r_j^2(x,\delta)$
但是上式的 Jacobian 矩阵很特殊：
$J(x,\delta) = \begin{bmatrix}\hat{J}&V\\0&D\end{bmatrix}$
对上面的式子用 Levenberg-marquardt method 方法求解即可。