Two-View 3D reconstruction using the sparse Levenberg–Marquardt algorithm
First, please consult this tutorial about the LM minimization and this document on the bundle adjustment.

Along with many other derivative based descent algorithms, the Levenberg-Marquardt algorithm relies on the partial derivative matrix, a.k.a. Jacobian Matrix, which is the matrix of all first-order partial derivatives of a vector-valued function:

$$
\begin{align}
J_i &amp;= \frac{\partial f(x_i, \beta)}{\beta} \\
\mathbf{J} &amp;= \begin{bmatrix}
    \dfrac{\partial f_1}{\partial x_1} &amp; \cdots &amp; \dfrac{\partial f_1}{\partial x_n}\\
    \vdots &amp; \ddots &amp; \vdots\\
    \dfrac{\partial f_m}{\partial x_1} &amp; \cdots &amp; \dfrac{\partial f_m}{\partial x_n} \end{bmatrix}
\end{align}
$$

Now how you compute Jacobian depends on what you actually want from the minimization. You could for example do a full bundle adjustment where Jacobian might have an involved expression. But, let me give an example of the simplest scheme I could think of.
Triangulation from multiple views is generally written in the form of minimization of the reprojection error. Here is the cost function for the two view problem:

$$
C(X) = d(x,\hat{x})^2 + d(x',\hat{x}')^2
$$
where $x$ and $x'$ are points in 2D views and $\hat{x}$ and $\hat{x}'$ are the projections of 3D points. In the matrix notation, the error function reads:

$$
\begin{pmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_m \\
\end{pmatrix} =
\begin{pmatrix}
\|x_1^a-P_a X^1\|+ \|x_1^b-P_b X^1\|\\
\vdots \\
\|x_m^a-P_a X^m\|+\|x_m^b-P_b X^m\| \\
\end{pmatrix}
$$
where $P_a, P_b$ are the projection matrices of cameras $a$ and $b$. Since we are searching for the 3D points only $\{\mathbf{X}_i\}$, the Jacobian is $Mx3$: 

$$
\begin{align}
\mathbf{J} &amp;= \begin{bmatrix}
    \dfrac{\partial \epsilon_1}{\partial x} &amp; \dfrac{\partial \epsilon_1}{\partial y} &amp; \dfrac{\partial \epsilon_1}{\partial z}\\
    \vdots &amp; \ddots &amp; \vdots\\
    \dfrac{\partial \epsilon_m}{\partial x} &amp; \dfrac{\partial \epsilon_m}{\partial y} &amp; \dfrac{\partial \epsilon_m}{\partial z} \end{bmatrix}
\end{align}
$$

Of course this problem requires a good initial solution, which could be obtained from a linear closed form solution. The quality of the non-linear minimizer is highly dependent on the initialization and the quality of 2D points.

For the rest, and different parameterizations, check out the links. 



Here are some remarks:


If you are not looking for re-inventing the wheel, then I would recommend Ceres Solver, which includes Automatic Differentiation. This saves the headache of Jacobians (and is much more accurate than numerical derivatives - numerical approximations are highly discouraged). Ceres also includes samples for bundle adjustment. 
If you are also optimizing for the pose (camera orientation), then I do not recommend operating on the 12 parameter vector of rotation and translations. Rather choose a representation such as quaternions or the exponential map (twist coordinates). As this preserves the manifold-ness, you will get much higher quality solutions, without ambiguities. Some examples are given in the links I provided.