Can we tackle the system directly?
Given $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}$ ($m$ equations in $n$ unknowns), how do we find a root?
Newton's method in multiple variables reads
$$x_{k+1} = x_k - Df(x_k)^{-1}f(x_k),$$with $Df : \mathbb{R}^{n} \rightarrow \mathbb{R}^{m\times n}$ is the Jacobian of $f$. It has quadratic convergence when starting close to a root when $Df(x^*)$ is invertible; linear convergence otherwise.
How do we invert the Jacobian at each step?
We are in fact approximating the solution via a Neumann series (for $x_0 = 0$, $\alpha = 1$)
$$x_n = \sum_{k=0}^{n-1}(I-A)^kb = P_{n-1}(A)b.$$Can we find a better polynomial approximation?
Ideally, we want a polynomial such that $P_{n}(\lambda_i) = \lambda_i^{-1}$ for all eigenvalues $\{\lambda_i\}_{i=1}^n$ of $A$.
Equivalently, we want to find $Q_n(A) = AP_{n}(A) - I$ which as roots at $\{\lambda_i\}_{i=1}^n$.
The Richardson iteration generates a solution in the Krylov subspace
$$K_n(A,b) = \text{span}\{b, Ab, A^2b, \ldots, A^{n-1}b\}.$$Can we find a better solution in this space?
Solve
$$\min_{x\in K_n(A,b)} \|Ax - b\|_2^2.$$N | Residual | Error | |
---|---|---|---|
20 | 1.80e-07 | 1.12e-06 | |
40 | 4.27e-07 | 5.53e-06 | |
60 | 7.67e-07 | 1.37e-05 | |
80 | 9.70e-07 | 2.20e-05 | |
100 | 9.53e-07 | 2.70e-05 |
What to do when we have more / less equations than unknowns?
The Moore-Penrose pseudo-inverse of a matrix $A$ is the unique matrix $A^\dagger$ satisfying:
Two special cases:
We can define the pseudo-inverse generally through the singular value decomposition:
$$A = U\Sigma V^*,$$$$A^\dagger = V_k \Sigma_k^{-1} U_k^*,$$with $k$ the rank of $A$.
How do we compute the minimum-norm solution of
$$\min_x \|Ax - b\|_2^2.$$The Jacobian satisfies
$$Df(\xi_k)(x_{k} - x_{k-1}) = f(x_k) - f(x_{k-1}),$$for $\xi_k$ an convex combination of $x_k, x_{k-1}$.
How do we get a usefull approximation $B_k$ or $H_k$ satisfying
$$H_k(x_{k} - x_{k-1}) = f(x_k) - f(x_{k-1}),$$or
$$(x_{k} - x_{k-1}) = B_k(f(x_k) - f(x_{k-1})),$$Assuming we have some $H_k$ (or $B_k$) satisfying the secant relation, how do we update it to obtain $H_{k+1}$ (or $B_{k+1}$)?
SR1-update:
$$H_{k+1} = H_k + \frac{(\Delta f_k - H_k\Delta x_k)(\Delta f_k - H_k\Delta x_k)^T}{(\Delta f_k - H_k\Delta x_k)^T\Delta x_k}$$$$B_{k+1} = B_k + \frac{(\Delta x_k - B_k\Delta f_k)(\Delta x_k - B_k\Delta f_k)^T}{(\Delta x_k - B_k\Delta f_k)^T\Delta f_k}$$with $g(x) = x - \alpha f(x)$.
which can be solved by
$$x_{k+1} = x_k - \alpha Df(x_k)^*\cdot f(x_k)$$Interdisciplinarity, own initiatives and originality are highly encouraged!