Trying to solve a set of non-linear equations

Question

I am looking for a solution to a set of non-linear equations. Let $K$ and $S$ be $N\times N$ real matrices. Given $K$ I would like to solve the following for $S$:

$$K_{jk} = \sum_a S_{ka} |S_{ja}|^\delta \mathrm{sign}(S_{ja}).$$

The exponent is $0 < \delta < 1$, but typically 1/2 and I can settle for a solution in that special case only, if that helps.

I realize that the solution will mostly likely be numerical. I have tried a brute-force solver in python, that converges well if $K$ is diagonal, but not really otherwise.

An example matrix where I'm having problems is $$K = \begin{pmatrix} 1.0 & 0.5 \\ 0.223 & 0.447 \end{pmatrix}.$$

I realize this is not a python forum, but here is my (ugly) code for reference:

import scipy.optimize
import numpy as np
N = 2
delta = 0.5
epsilon = 1e-3
K = np.array([[1.0, 0.5],
              [0.223, 0.447]])
x0 = K.flatten()**0.666
def f(x):
    return np.array([[sum([x[k*N+a]np.abs(x[jN+a])*deltanp.sign(x[j*N+a]) for a in range(N)])-K[j, k] for k in range(N)] for j in range(N)]).flatten()
def fprime(x):
    ret = np.array([[np.array([[sum([(kdN+jd == kN+a)np.abs(x[jN+a])deltanp.sign(x[jN+a]) + ((kdN+jd == jN+a) and (np.abs(x[j*N+a]) > epsilon))x[kN+a]deltanp.abs(x[j*N+a])(delta-1.)np.sign(x[jN+a]) for a in range(N)])-K[j, k] for k in range(N)] for j in range(N)]) for jd in range(N)] for kd in range(N)])
    return ret.reshape(len(x), -1)
S, infodict, status, msg = scipy.optimize.fsolve(f, x0, fprime=fprime, xtol=1e-16, epsfcn=1e-16, maxfev=2000, full_output=True)
printmsg = f"{msg} nfev = {infodict['nfev']}, njev = {infodict['njev']}"
assert status == 1, printmsg

Any feedback is welcome.

PierreCarre · Answer 1 · 2021-08-18T10:01:48.333

What was exactly the brute force method you used? I tried on Wolfram Mathematica (Findroot[]) with $$ K = \left( \begin{array}{ccc} 4 & 1 & 0 \\ 1 & 4 & 1 \\ 0 & 1 & 4 \\ \end{array} \right) $$

and got

$$ S = \left( \begin{array}{ccc} 2.50426 & 0.111121 & -0.000214794 \\ 0.10846 & 2.48975 & 0.10846 \\ -0.000214794 & 0.111121 & 2.50426 \\ \end{array} \right) $$ with a residual within machine precision.

Not very pretty Mathematica code:

(* Defining the matrix  K *)
K = Table[
   If[i == j, 4, If[Abs[i - j] == 1, (-1)^(i + j), 0]], {i, 1, 3}, {j,
     1, 3}];
(* K={{1.,0.5},{0.223,0.447}};*)
n = Length[K];
(* Solving the nonlinear system *)
num = FindRoot[
   Table[
    K[[j, k]] - 
     Sum[S[k, a]Abs[S[j, a]]^(1/2)Sign[S[j, a]], {a, 1, n}], {j, 1, 
     n}, {k, 1, n}], 
   Flatten[Table[{S[j, k], K[[j, k]]}, {j, 1, n}, {k, 1, n}], 1], 
   Method -> {"Newton", "StepControl" -> "TrustRegion"}];
(* Building the solution matrix S *)
sol = Table[S[i, j] /. num, {i, 1, n}, {j, 1, n}];
(* Displaying the solution *)
sol // MatrixForm

Thanks Pierre, I replied by updating the problem description. My matrix is very simple, but my algo must be very bad. I basically break this down to an $N\times N$ variable root finding problem. If you can share your Mathematica code (I don't know the first thing about the language) then I can try myself. Actually if it could also run on Wolfram Alpha that would be easier... — poldpold, Aug 18 '21 at 09:19

Trying to solve a set of non-linear equations

1 Answers1