2
$\begingroup$

In a personal project of mine, I've derived a couple of equations.

$$ \sum_{j=1}^Nu_{ij}^*(\theta_1,...,\theta_N,J) - k_{i} = 0 \qquad \forall 1\leq i \leq N$$ $$\sum_{i<j}^N\big(u_{ij}^*(\theta_1,...,\theta_N,J)\big)^2 - O^* = 0$$

where $u_{ij}^*(\theta_1,...,\theta_N,J)$ is a solution to the equation

$$u_{ij} =\frac{1}{2} + \frac{1}{2}\tanh{\left( - \frac{\theta_i + \theta_j}{2} + 2\cdot J\cdot M \cdot u_{ij}\right)} $$

where $M,N,k_i^*, O^*$ are known parameters, and $\theta_1,...,\theta_N,J$ are $N+1$ unknown parameters. I thought this system should be solvable since the top two equations give us $N+1$ equations as well.

I think the difficulty in this problem is due to the presence of some unknown/hidden function $u_{ij}^*$.

My first try was to construct a quadratic cost function, and to minimize this quadratic cost function using gradient descent (involving the partial derivatives of $u_{ij}^*$ (which can be obtained implicitly), however this method did not work out.

Does anyone have any idea on how to solve this or know of some tools which I can use to solve this?

Or maybe you could point me to a different place?

$\endgroup$
  • 1
    $\begingroup$ What do you mean "this method did not work out"? Was the method inadequate, or was there a bug in your code? $\endgroup$ – Wolfgang Bangerth Jun 20 '17 at 15:28
  • $\begingroup$ I think that you have more equations than unknowns. Your unknowns are $\mathbf{x} = (u_{1,1}, u_{1,2}, \cdots, u_{i,j},\cdots, u_{N-1, N}, u_{N,N}, \theta_{1}, \theta_{2}, \cdots, \theta_{N})$, that add up to $N^2 + N$. While you have $N^2 + N + 1$ equations. $\endgroup$ – nicoguaro Jun 20 '17 at 16:48
  • $\begingroup$ I solved using scipy.optimize.root (neglecting the equation that involves the $O*$) and using your approach with scipy.optimize.minimize for $N=2$, and it worked. Of course, some of the components of the solution vector are different. $\endgroup$ – nicoguaro Jun 20 '17 at 17:05
  • $\begingroup$ @WolfgangBangerth , depending on my initial parameter value choice ( in the case of convergence) it converged to local minima of the cost function. $\endgroup$ – hoolaboris Jun 21 '17 at 10:03
  • $\begingroup$ I see, I missed $J$. I will post an answer considering that. $\endgroup$ – nicoguaro Jun 21 '17 at 17:39
2
$\begingroup$

First, let's rewrite your system in the form of a vector function for which you want to find a root.

$$\mathbf{F}(\mathbf{x}) = 0$$

with

\begin{align} F_i(\mathbf{x}) &= \sum_{j=1}^N u_{ij} - 1 \qquad \forall 1\leq i \leq N\\ F_i(\mathbf{x}) &= \sum_{i<j}^N u_{ij}^2 k_i^2 - R \qquad i = N + 1\\ F_i(\mathbf{x}) &= u_{ij} -\frac{1}{2} - \frac{1}{2}\tanh{\left( - \frac{\theta_i + \theta_j}{2} + 2 J u_{ij}\right)} \qquad \forall i >N +1 \land i\leq N^2 + N + 1\, , \end{align}

and

$$\mathbf{x} = (u_{1,1},u_{1,2},\cdots,u_{i,j},\cdots, u_{N-1,N},u_{N,N},θ_1,θ_2,\cdots,θ_N,J)\, .$$

The parameters $k_i$ are absorbed as part of the new variables $u_{ij}$, and the same happened to $M$. That way, you have $N + 1$ less parameters to worry about (I also renamed $O^*$ as $R$).

Then, you can use a root finding algorithm. If you can compute the Jacobian, you can use something like Newton's method.

The option that you propose would be something like

$$\min_{\mathbf{x}} \Vert \mathbf{F}(\mathbf{x})\Vert^2\, ,$$

what you can plug into a minimization algorithm

Below, I implemented both approaches in Python using scipy.optimize. Notice that basin-hopping is a global minimization approach and it is not a fair comparison. Also, note that basin-hopping will run several initial points and it will take more time that just using a single initial point. If you have a good guess, then it might improve.

from __future__ import division, print_function
import numpy as np
from numpy.linalg import norm
from scipy.optimize import root, basinhopping


def fun(x, R, N, k):
    u = x[:N**2].copy()
    u.shape = (N, N)
    f = np.zeros_like(x)
    f[:N] = -np.ones((N))
    f[N] = -R
    for row in range(N):
        for col in range(N):
            f[row] += u[row, col]
            f[N + 1 + N*col + row] = u[row, col] - 0.5 -\
                0.5*np.tanh(-0.5*(x[N**2 + row] + x[N**2 + col]) +
                2*x[-1]*u[row, col])
            if row < col:
                f[N] += u[row, col]**2*k[row]**2
    return f


def opt_fun(x, R, N, k):
    return norm(fun(x, R, N, k))**2


N = 2
R = 1
size = N**2 + N + 1
x0 = np.random.normal(loc=1,size=size)
k = np.ones((N))
res_root = root(fun, x0, args=(R, N, k))
res_min = basinhopping(opt_fun, x0,
                       minimizer_kwargs={'args': (R, N, k)}, seed=1)
print("Root finding algorithm norm: {:.6g}".format(norm(res_root.fun)))
print("Minimization algorithm norm: {:.6g}".format(norm(res_min.fun)))

For $N=2$ the root finding works better. For higher orders it seems that the optimization procedure works better. For example for $N=10$, this is the output

Root finding algorithm solution norm: 0.772523
Minimization algorithm solution norm: 1.67062e-13
$\endgroup$
  • $\begingroup$ Thanks for the answer. There are some mistakes in your code where you fill in the vector function, but I understand the general idea. I'm however not sure about your redefinition of $u_{ij}$, as it does not seem consistent between the first $N$ equations and the $N+1$'th equation $\endgroup$ – hoolaboris Jun 22 '17 at 12:04
  • $\begingroup$ What are the mistakes? And what is the inconsistency in $u_{ij}$? If you point them out I can correct them or discuss them with you. $\endgroup$ – nicoguaro Jun 22 '17 at 12:39
  • $\begingroup$ According to your equations posted above the code: in the code you're adding -1 to the first $N$ equations for every column, i think it should only be added before the 2nd loop begins. For the $N+1$'th equation, you are again adding $R$ for every column/row, while it should only be added once. For the equation $i$ where $(N+1)< i \leq N^2+N+1$, the index is wrong. I think it should be $N+1 + N*col + row$ instead of $N+1+col*row + row$. $\endgroup$ – hoolaboris Jun 22 '17 at 13:33
  • $\begingroup$ For the inconsistency on $u_{ij}$, i'll get back to you later :) But i think the redefinition is fine in the first $N$ equations, but in the $(N+1)$'th equation you can't use the same redefinition for $u_{ij}$, since you now sum over $i$ aswell (which kind of changes the system of equations), but I'll think about it some more later and get back to you! $\endgroup$ – hoolaboris Jun 22 '17 at 13:39
  • $\begingroup$ @hoolaboris, I think that you are right. I have fixed it. I also changed the initial guess for a random one. $\endgroup$ – nicoguaro Jun 22 '17 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.