How to substitute numeric values in a symbolic Jacobian matrix?

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.

-
What exactly do you intend J[X] = D[f[X], {X}]; to do? – Mr.Wizard Jun 14 '12 at 17:19
MJ asking, JM editing ... Hmmmm – Dr. belisarius Jun 14 '12 at 17:30
@belisarius Don't you know that it's write your name in reverse day? ‮ – Mr.Wizard Jun 14 '12 at 17:45
By J[X] we are trying to find the Jacobian matrix. This is the way in defining Newton root-finding method for solving nonlinear systems. Please also see the first example in the "Application" part of "LinearSolve" tag in the help of Mathematica (documentation center). – M.J. Jun 14 '12 at 18:08
M.J. let me phrase that differenty: you seem to be attempting to create a function but you are not, therefore J[Y0] does not evaluate. – Mr.Wizard Jun 14 '12 at 18:25

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}]]},
]

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

-
Another way would be to use indexed variables, e.g. D[{Sin[C[1]] - Cos[C[2]], C[2] Cos[C[1]]}, {{C[1], C[2]}}] generates an appropriate Jacobian. – J. M. Jun 14 '12 at 23:03
@J.M. I did not know that variables for D need not be a symbol. This is good to know. – Michael Wijaya Jun 14 '12 at 23:07
The descriptions are so much helpful. Thanks. – M.J. Jun 15 '12 at 5:21
There is still a big problem. For large dim, e.g. dim=10000, the computation of the Jacobian in this way takes so much time! I used SparseArray[] to reduce the time, but coud not succeed. Please give some tips for this action. – M.J. Jun 15 '12 at 7:19
There is a still a problem in this way, you obtain the Jacobian at one point, what should we do when we wish to call it regularly in a loop? – M.J. Jun 15 '12 at 10:40

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_] :=


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.

-
Thanks a lot Jen, it works excellent now. – M.J. Jun 15 '12 at 5:21
There is still a big problem. For large dim, e.g. dim=10000, the computation of the Jacobian in this way takes so much time! I used SparseArray[] to reduce the time, but coud not succeed. Please give some tips for this action. – M.J. Jun 15 '12 at 6:53
I don't see a big problem: just replace the definition of jacobiMatrix by this: jacobiMatrix = D[SparseArray[f[X]], {X}]; – Jens Jun 15 '12 at 7:26
Please do it, in this way and based on applying SparseArray[] on f[X], the Jacobian will not take any numeric values to itself! Please Test it. Consider that we are next to solve LinearSolve[J[Y0], f[Y0]], the system says that there is no solution! – M.J. Jun 15 '12 at 7:38
I think SparseArray[] does not work properly with the command /.Thread... – M.J. Jun 15 '12 at 7:39