How to substitute numeric values in a symbolic Jacobian matrix?

Question

I have a multi-variate function from $\mathbb{R}^n\to\mathbb{R}^n$. Choosing any desired initial vector, we can produce the corresponding function value, which is a vector as follows. The main problem lies in the computation of the Jacobian matrix. In the code below:

dim = 10;
X = Table[Subscript[x, i], {i, 1, dim}]
f[x_] := (n = Length[x]; 
   Table[If[i <= n - 1, x[[i]] x[[i + 1]] - 1, x[[1]] x[[n]] - 1], {i,
      1, n}]);
f[X]
J[X] = D[f[X], {X}];
J[X]
Y0 = Table[2, {i, 1, dim}];
f[Y0]
J[Y0]

the Jacobian (the differentiation of f[x]) is defined (theoretically) correctly, but I cannot plug any data to obtain the numerical values of the Jacobian matrix!? It would be appreciated if anyone take a look and give some suggestions to find J[Y0] correctly? Here Y0 is the initial vector chosen arbitrary.

Answer 1

From what I understand, you want to compute the Jacobian matrix of a function $f$ at $(x_1,x_2,\ldots,x_n)$ and evaluate it at a given point, say $(2,2,\dots,2)$.

This seems to be the function you are interested in:

f[{X__}] := {X} RotateLeft[{X}, 1] - 1;

In[2]:= f[{x1,x2,x3}]
Out[2]= {-1+x1 x2,-1+x2 x3,-1+x1 x3}

Since the dimension may change, we want to use SetDelay or := so that the Jacobian function recomputes the matrix every time it is called.

J[{X__}] := D[f[{X}], {{X}}]

In[4]:= J[{x1,x2,x3}]
Out[4]= {{x2,x1,0},{0,x3,x2},{x3,0,x1}}

If you want the Jacobian matrix of $f$ at the point $(2,3,4)$, you can do the substitution

In[5]:= J[{x1,x2,x3}] /. Thread[{x1,x2,x3}->{2,3,4}]
Out[5]= {{3,2,0},{0,4,3},{4,0,2}}

I would also recommend that you avoid using subscripts. Here is one way to generate a vector of variables such as {x1,x2,x3,x4,x5} programatically:

In[6]:= ToExpression[Table["x"~~ToString[k],{k,1,5}]]
Out[6]= {x1,x2,x3,x4,x5}

variableVector[n_] := ToExpression[Table["x" ~~ ToString[k], {k, 1, n}]];

Now we can write a function which computes the Jacobian at a given point as follows:

Jvalue[{Y__}] := With[
  {X=variableVector[Length[{Y}]]},
  J[X] /. Thread[X->{Y}]
 ]

In[9]:= Jvalue[{2,3,4}]
Out[9]= {{3,2,0},{0,4,3},{4,0,2}}

Update

Here is a a function which produces the Jacobian matrix at a given vector:

sparse[v_] := (
  Clear[x];
  xVariables = Map[x, Range[10000]];
  Evaluate[xVariables] = v;
  SparseArray[{{10000, 10000} -> x[1], {10000, 1} -> x[10000], {i_, i_} -> x[i + 1], ({i_, j_} /; j - i == 1) -> x[i]}, {10000, 10000}]
 )

Now you should be able to use it in a loop:

In[2]:= sparseTest1=sparse[RandomInteger[{1,20},10000]];
        sparseTest2=sparse[RandomInteger[{1,20},10000]];

In[3]:= sparseTest1[[1,2]]
Out[3]= 11

In[4]:= sparseTest2[[1,2]]
Out[4]= 8

Answer 2

To pinpoint exactly where your approach fails, I'll stay as close as possible to what you did. The two lines I changed are separated from your code by a blank line.

dim = 10;
X = Table[Subscript[x, i], {i, 1, dim}]
f[x_] := (n = Length[x]; 
   Table[If[i <= n - 1, x[[i]] x[[i + 1]] - 1, x[[1]] x[[n]] - 1], {i,
      1, n}]);
f[X]

jacobiMatrix = D[f[X], {X}]
J[xVariable_] := jacobiMatrix /. Thread[X -> xVariable]

Y0 = Table[2, {i, 1, dim}];
f[Y0]
J[Y0]

All I did is to separate the definition of the symbolic Jacobi matrix jacobiMatrix from the definition of the function that you want to call J[X].

In the function definition, I purposely named the vector argument xVariable instead of X, to make it clear that this is a "dummy" argument that can take numerical (or other) values, and not the same as the symbolic vector X you already defined.

What the function does on the right-hand side is to substitute the instantaneous values passed to it via xVariable in place of the symbolic entries $x_1, x_2\ldots$ contained in J. The Thread is used to make the rule arrow -> apply individually to each pair of element in X and xVariable.

In your original definition of J[X], you didn't actually define a function because X already had been assigned a value, which was inserted at the time J[X] was defined. If you use your definition and then inspect what has been defined by doing ?J, you'll see that Mathematica knows the value of J[X] as you stated it, but only if X is literally equal to $\{x_1, x_2\ldots\}$. The crucial thing if you want to use numerical variables is to replace these formal variables in the formal matrix. That's what the line J[xVariable_] in my modified code does.

Edit

If the Jacobi matrix is so big that you need to use SparseArray to represent it, the rules in the Thread command can be applied to the ArrayRules instead. The only change is this:

jacobiMatrix = D[SparseArray[f[X]], {X}]
J[xVariable_] := 
 SparseArray[ArrayRules[jacobiMatrix] /. Thread[X -> xVariable]]

Now when you calculate J[Y0], the result will be a SparseArray with the correct entries. To check the entries, you could say Normal[J[Y0]], but for a large matrix that's of course a crazy idea. Instead, one could do

ArrayPlot[J[Y0]]

to visualize the result.