# The most optimal way to solve this set of non-linear equations in high dimensions

So I have a series of non-linear equations which I wish to solve as fast as possible, to illustrate for the case of $n = 4$, I have the following equations:

\begin{gather*} c_0\exp(\lambda_0-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = 1 \\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_1 \\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_2 \\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_3 \\ c_4\exp(\lambda_0+\lambda_4-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_4 \end{gather*} where \begin{gather*} f = \int_{d_4}^{\infty} \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_0 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_1 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_2 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_3 = \int_{-\infty}^{d_4} \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_4 = \int_{d_4}^{\infty} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ \end{gather*}

The definitions are as follows:

• $P_1, P_2, P_3, P_4$ are known values.

• $d_1, d_2, d_3, d_4$ are known values.

• $q(x_1, x_2, x_3, x_4)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.

• Thus $f, c_0, c_1, c_2, c_3, c_4$ are all simply real numbers.

For the general case I have:

\begin{gather*} c_0\exp(\lambda_0-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = 1 \\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_1 \\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_2 \\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_3 \\ \vdots \\ c_n\exp(\lambda_0+\lambda_n-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_n \end{gather*} where \begin{gather*} f = \int_{d_n}^{\infty} \cdots \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_0 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_1 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_2 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_3 = \int_{-\infty}^{d_n} \cdots \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ \vdots \\ c_n = \int_{d_n}^{\infty} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ \end{gather*}

The definitions are as follows:

• $P_1, P_2, P_3, \cdots, P_n$ are known values.

• $d_1, d_2, d_3, \cdots, d_n$ are known values.

• $q(x_1, x_2, x_3, \cdots, x_n)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.

• Thus $f, c_0, c_1, c_2, c_3, \cdots, c_n$ are all simply real numbers.

At most $n$ will be 35, so there will be 36 equations to solve simultaneously. What is the best method to solve this as quickly as possible using whatever programming language possible? Also I will need to change parameter inputs for $P_i$, $d_i$ and the variance covariance matrix and solve for a new set of $\lambda$'s, I will need to do this at least 2000 times, which means effectively I will need to solve a set of 35 equations for 2000 times.

The first is Newton's method with a line-search. Basically, you'll have a system of equations $F(x)=0$ and you'll be solving iterates $F^\prime(x) \delta x = -F(x)$ where $F^\prime(x)$ is the Frechet (total) derivative of $F$. It's a linear operator, so you can solve this in something like Matlab or use LAPACK in C++ or one of the numerical libraries in Python. Then, you'll have to do a line-search based on some merit function like $\frac{1}{2}\|F(x+\alpha \delta x)\|^2$ for $\alpha$ and take a step $x+\alpha \delta x$. The line-search is necessary in order to guarantee that Newton's Method converges.
If you have an equality constrained optimization code lying around, you can cheat this and solve $$\min_{x,y} \{ \frac{1}{2} \|y\|^2 : F(x) = y \}$$ Any good optimization code will do the globalization for you, so it's normally easiest to go this route assuming you have a constrained optimization code.