[NLsolve] Providing the Jacobian of a system when variables are provided in an Array?

New to Julia
1

When solving systems of equation, sometimes, it’s better to have the variables in Array format, instead of a vector.

That said, can we provide multidimensional Jacobians to the solvers of the NLsolve package?

For example, I want to solve the system:

\begin{array}{r} 3 x_{11}^{3}+2 x_{11} x_{21}=5 \\ x_{21}^{3}+2 x_{22}=3 \\ 2 x_{12}+3 x_{22}^{3}=5 \\ 2 x_{12}+3 x_{21}=5 \end{array}

The variables are represented in array format (instead of a vector).

x=\left[\begin{array}{ll} x_{11} & x_{12} \\ x_{21} & x_{22} \end{array}\right]

The solution is:

x^{*}=\left[\begin{array}{ll} 1 & 1 \\ 1 & 1 \end{array}\right]

Solving the system in vector format looks like this.

using NLsolve

# solution [1,1,1,1]

function f!(F, x)
    F[1]= 3*x[1]^3  + 2*x[1]*x[2]   - 5.
    F[2]= x[2]^3    + 2*x[4]        - 3.
    F[3]= 2*x[3]    + 3*x[4]^3      - 5.
    F[4]= 2*x[3]    + 3*x[2]        - 5.
end

function j!(J, x)
    J[1, 1] = 9*x[1]^2 + 2*x[2]
    J[1, 2] = 2*x[1]
    J[1, 3] = 0.
    J[1, 4] = 0.

    J[2, 1] = 0.
    J[2, 2] = 3*x[2]^2
    J[2, 3] = 0.
    J[2, 4] = 2.

    J[3, 1] = 0.
    J[3, 2] = 0.
    J[3, 3] = 2.
    J[3, 4] = 9*x[4]^2

    J[4, 1] = 0.
    J[4, 2] = 3.
    J[4, 3] = 2.
    J[4, 4] = 0.
end

nlsolve(f!,j!,[0., 0., 0., 0.])

Solving it while representing the variables in an array looks like this.

## system in "matrix form"

function f!(F, x)
    F[1,1] = 3*x[1,1]^3  + 2*x[1,1]*x[2,1]   - 5.
    F[2,1] = x[2,1]^3    + 2*x[2,2]          - 3.
    F[1,2] = 2*x[1,2]    + 3*x[2,2]^3        - 5.
    F[2,2] = 2*x[1,2]    + 3*x[2,1]          - 5.
end

nlsolve(f!, [0. 0.; 0. 0.])

Is there any way to provide the Jacobian to NLsolve in the second case?

2

Why don’t you pass a vector of parameters and reshape it to a matrix within f!? This way the Jacobian is still a matrix. You may reshape the vector to a matrix once you have your solution

3

For convenience. The actual inputs of my non-linear system are multidimensional arrays that represent the value functions of different states of the world. To reshape that is very cumbersome.

4

If you “flatten” the vectors with vec, and then reshape with reshape, it should yield the same results. Although, NLsolve takes any AbstractArray.

5

@amrods, yes… NLsolve takes any AbstractArray and works fine with arrays as inputs if you don’t provide the Jacobian. I didn’t see anything in the Docs about multidimensional Jacobians tho.

6

What if you supply the Jacobian with autodiff? Even if you really want to supply analytical expressions, it can serve you to compare what is the structure of a multidimensional Jacobian.

7

The structure with AD is a matrix for the Jaobian and veced input. I’m not sure there’s a another/better way to do it, because the linear algebra has to work out. Currently, all input is veced where it’s used. If it’s an AbstractVector it’s a noop. There no similar mat that makes it NxM as far as I’m aware, and it’d be ambiguous either way because there are potentially many possible matrices that fit the number of elements.

The “input a containter of a different shape as your variables” is a convenience layer. It doesn’t work well for Jacobians. We have the same in Optim where you can input non-Vector variables and gradients, but if you want to use a Hessian-based optimizer, you have to use a Matrix.

2 Likes