9

Problem

Solving a non-linear system of equations.

The number of variables is the same as the number of equations.

When I fix a set of variables (say $\vec{y}$) and keep another set free (say $\vec{x}$), the system becomes an under-determined, dense, and linear system of the subset of variables $$A(\vec{y})\vec{x} = \vec{b}(\vec{y})\tag{1}\label{system1}$$

where $A(\vec{y})$ is a dense matrix, and $\vec{b}(\vec{y})$ is a dense vector). Let's call this sub-system $\eqref{system1}$ as system 1.

When I fix $\vec{x}$ and keep $\vec{y}$ free, the system becomes an over-determined, sparse, and non-linear system of the subset of variables $$F(\vec{x}, \vec{y}) = 0\tag{2}\label{system2}$$

The Jacobian $J(\vec{x}, \vec{y})$ is available in a closed form. Let's call this sub-system $\eqref{system2}$ as system 2.

About half of the equations in $\eqref{system2}$ are equality constraints that are linear in terms of $\vec{y}$. Each of the constraints is quite sparse and involves only about 5% of all variables.

Can I solve it with the following?

Algorithm 1

  1. Initialize $\vec{x} = \vec{x}_0$, and $\vec{y} = \vec{y}_0$
  2. Fix $\vec{y}_{n - 1}$, solve $A(\vec{y}_{n-1})\vec{x}_{n} = \vec{b}(\vec{y}_{n-1})$ for one of all the possible $\vec{x}_{n}$ because this system is under-determined.
  3. Fix $\vec{x}_{n - 1}$. Perform one iteration of Newton's method for solving $F(\vec{x}_{n-1}, \vec{y}_{n}) = 0$ for $\vec{y}_{n}$.
  4. If not converged, go to step 2.

Algorithm 2

If I replace step 2 in Algorithm 1 by a Newton's method iteration for solving system 1 $\eqref{system1}$, then I guess the steps become a block Newton's method.

Question

But I don't know if these two algorithms can work because system 1 $\eqref{system1}$ is under-determined and system 2 $\eqref{system2}$ is over-determined.

Can this work?

Related

R zu
  • 163
  • 9
  • How do you choose which variables belong to set $\vec{x}$ and which to $\vec{y}$? – Yuriy S Oct 10 '18 at 15:56
  • The system is linear in terms of the variables $\vec{x}$. The remaining variables are $\vec{y}$. – R zu Oct 10 '18 at 16:04
  • I see. What I would do at a first glance, is add however many variables from Y to the X subset to make the number of equations match number of variables, and linearize the corresponding equations to get the first approximation. Then use the Newton's method on both subsets as intended. But that's just an idea, I'm not sure if it will work – Yuriy S Oct 10 '18 at 16:07
  • Can't linearize (theoretically without numerical methods) around the solution because I don't have a good theoretical approximation of the solution. I think Newton's method is equivalent to linearization at the current approximation of the solution, using the Jacobian. Moving variables from Y to X has several disadvantages in terms of programming. It makes system 1 not linear and not dense anymore. That means higher memory usage and forced to use slower methods. Solving sparse system is much slower than solving a dense system. Speed is crucial as I have to speed up the program by 50 times. – R zu Oct 10 '18 at 16:20
  • @YuriyS Maybe. Suppose we linearize the whole thing around the current approximation of the solution. Then the question is equivalent to solving a linear system by block matrix method with the variables split unevenly. If that can work, then the scheme in the question can work. – R zu Oct 10 '18 at 16:44
  • Use automatic differentiation to linearize? – Chris Rackauckas Oct 15 '18 at 14:52
  • @ChrisRackauckas I already have the Jacobian in close form (I do it by hand, no need for automatic differentiation). Newton's method is equivalent to linearization, using the Jacobian. Also, linearize around what point? How would linearization prove/disprove that the algorithm would work? – R zu Oct 15 '18 at 15:06
  • It's similar to iterative least squares approximations in ML. It can be hard to prove that it will converge, but even so it's useful in practice. I suspect similarly for this: no one can prove/disprove that it will work, but with a decent initial condition you should be fine. It's worth a shot at least. – Chris Rackauckas Oct 15 '18 at 15:15
  • Can you tell me a bit more about the iterative least squares that you refer to? When I search online for iterative least squares, there are too many different things. For convergence, some algorithm needs positive definite matrix to work, and some depends on the eigenvalues of the Jacobian. I don't know how well would algorithm 1/2 work. – R zu Oct 15 '18 at 15:18
  • As your bounty is almost over, maybe you can think of what kind of answer you really want? Because I doubt anyone can answer "can this work" without first trying/doing some experiments with the method, using your data. At least that would be the most useful way to check the validity of the method. I'm sure there's some general theory available, but checking directly is usually my first choice – Yuriy S Oct 16 '18 at 13:36
  • I thought it is "just" linear algebra. If no one can see obvious loopholes of this method (about convergence and stability), I guess I will just go ahead and program the whole thing. I might also post the question in sundial/other non-linear optimization/solver's forum. Or try theoretical computer science stack exchange. – R zu Oct 16 '18 at 14:57
  • Probably can work because this is equivalent to a block triangular approximation of the Jacobian. I got the answer from the SUNDIALS mail list: http://sundials.2283335.n4.nabble.com/Solving-a-non-linear-system-with-a-linear-sub-system-td4654414.html – R zu Oct 19 '18 at 19:39
  • Probably can work. This is equivalent to a block triangular approximation of the Jacobian. I am still not completely sure if the unbalance of variables would cause any problem. – R zu Oct 19 '18 at 19:39

0 Answers0