diff --git a/nlls/lm_test.go b/nlls/lm_test.go
new file mode 100644
index 0000000..4dbf22f
--- /dev/null
+++ b/nlls/lm_test.go
@@ -0,0 +1,245 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package nlls
+
+import (
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/mat"
+	"gonum.org/v1/gonum/stat/distuv"
+)
+
+// LeastSquares is a type for solving linear least squares problems with LM.
+type LeastSquares struct {
+	X *mat.Dense
+	Y []float64
+}
+
+func (l LeastSquares) Problem() LMProblem {
+	r, c := l.X.Dims()
+	return LMProblem{
+		Dim:  c,
+		Size: r,
+		Func: l.Func,
+		Jac:  l.Jac,
+	}
+}
+
+func (l LeastSquares) Optimal() []float64 {
+	r, c := l.X.Dims()
+	if len(l.Y) != r {
+		panic("size mismatch")
+	}
+	yVec := mat.NewVecDense(len(l.Y), l.Y)
+
+	tmp := make([]float64, c)
+	tmpVec := mat.NewVecDense(len(tmp), tmp)
+
+	err := tmpVec.SolveVec(l.X, yVec)
+	if err != nil {
+		panic("singular")
+	}
+	return tmp
+}
+
+func (l LeastSquares) Func(dst, params []float64) {
+	l.funcJac(nil, dst, params)
+}
+
+func (l LeastSquares) Jac(dst *mat.Dense, params []float64) {
+	l.funcJac(dst, nil, params)
+}
+
+func (l LeastSquares) funcJac(jacDst *mat.Dense, funDst, params []float64) {
+	if funDst != nil {
+		for i := 0; i < len(funDst); i++ {
+			x := l.X.RawRowView(i)
+			diff := floats.Dot(x, params) - l.Y[i]
+			funDst[i] = diff
+		}
+	}
+	if jacDst != nil {
+		jacDst.Copy(l.X)
+	}
+}
+
+func constructLeastSquares(trueParam []float64, noise float64, offset bool, nData int, source rand.Source) *LeastSquares {
+	norm := rand.New(source).NormFloat64
+	dim := len(trueParam)
+	xs := mat.NewDense(nData, len(trueParam), nil)
+	ys := make([]float64, nData)
+	for i := 0; i < nData; i++ {
+		if offset {
+			xs.Set(i, 0, 1)
+		} else {
+			xs.Set(i, 0, norm())
+		}
+		for j := 1; j < dim; j++ {
+			xs.Set(i, j, norm())
+		}
+
+		x := xs.RawRowView(i)
+		y := floats.Dot(trueParam, x) + distuv.Normal{Mu: 0, Sigma: noise, Src: source}.Rand()
+		ys[i] = y
+	}
+	return &LeastSquares{
+		X: xs,
+		Y: ys,
+	}
+}
+
+// Powell, M. J. D. "A Hybrid Method for Nonlinear Equations", in P. Rabinowitz, ed.,
+// "Numerical Methods for Nonlinear Algebraic Equations", Gordon and Breach, 1970.
+func powellFunc(dst, x []float64) {
+	dst[0] = x[0]
+	dst[1] = 10*x[0]/(x[0]+0.1) + 2*x[1]*x[1]
+}
+
+func powellJac(dst *mat.Dense, x []float64) {
+	dst.Set(0, 0, 1)
+	dst.Set(0, 1, 0)
+	dst.Set(1, 0, math.Pow(x[0]+0.1, -2))
+	dst.Set(1, 1, 4*x[1])
+}
+
+// The following test functions are taken form:
+// - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained optimization software. 
+//   ACM Trans Math Softw 7 (1981), 17-41.
+func bealeFunc(dst, x []float64) {
+	dst[0] = 1.5 - x[0]*(1-x[1])
+	dst[1] = 2.25 - x[0]*(1-math.Pow(x[1], 2))
+	dst[2] = 2.625 - x[0]*(1-math.Pow(x[1], 3))
+}
+
+func biggsEXP6Func(dst, x []float64) {
+	for i := 0; i < 13; i++ {
+		z := float64(i) / 10
+		y := math.Exp(-z) - 5*math.Exp(-10*z) + 3*math.Exp(-4*z)
+		dst[i] = x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) + x[5]*math.Exp(-x[4]*z) - y
+	}
+}
+
+func extendedRosenbrockFunc(dst, x []float64) {
+	dst[0] = 10*(x[1]-x[0]*x[0])
+	dst[1] = 1-x[0]
+}
+
+type LMTest struct {
+	prob     LMProblem
+	expected []float64
+	tol      float64
+	settings Settings
+}
+
+func TestLM(t *testing.T) {
+	// Constract a linear least squares probelm.
+	trueParam := []float64{0.7, 0.8, 5.6, -40.8}
+	nData := 50
+	noise := 1e-2
+	ls := constructLeastSquares(trueParam, noise, true, nData, rand.NewSource(1))
+	optLSParams := ls.Optimal()
+
+	// Numerical Jacobians for problems.
+	bealeNumJac := NumJac{Func: bealeFunc}
+	biggsNumJac := NumJac{Func: biggsEXP6Func}
+	rosenbrockNumJac := NumJac{Func: extendedRosenbrockFunc}
+
+	problems := []LMTest{
+		// Simple linear fit problem.
+		LMTest{
+			prob: LMProblem{
+				Dim:        len(trueParam),
+				Size:       nData,
+				Func:       ls.Func,
+				Jac:        ls.Jac,
+				InitParams: nil,
+				Tau:        1e-3,
+				Eps1:       1e-8,
+				Eps2:       1e-8,
+			},
+			expected: optLSParams,
+			tol:      1e-6,
+			settings: Settings{Iterations: 100, ObjectiveTol: 1e-16},
+		},
+		// Powell problem.
+		LMTest{
+			prob: LMProblem{
+				Dim:        2,
+				Size:       2,
+				Func:       powellFunc,
+				Jac:        powellJac,
+				InitParams: []float64{3, 1},
+				Tau:        1,
+				Eps1:       1e-15,
+				Eps2:       1e-15,
+			},
+			expected: []float64{0.0, 0.0},
+			tol:      1e-2,
+			settings: Settings{Iterations: 100, ObjectiveTol: 1e-16},
+		},
+		// Beale problem.
+		LMTest{
+			prob: LMProblem{
+				Dim:        2,
+				Size:       3,
+				Func:       bealeFunc,
+				Jac:        bealeNumJac.Jac,
+				InitParams: []float64{1, 1},
+				Tau:        1,
+				Eps1:       1e-15,
+				Eps2:       1e-15,
+			},
+			expected: []float64{3.0, 0.5},
+			tol:      1e-5,
+			settings: Settings{Iterations: 100, ObjectiveTol: 1e-16},
+		},
+		// Biggs EXP6 problem.
+		LMTest{
+			prob: LMProblem{
+				Dim:        6,
+				Size:       13,
+				Func:       biggsEXP6Func,
+				Jac:        biggsNumJac.Jac,
+				InitParams: []float64{1, 2, 1, 1, 1, 1},
+				Tau:        1e-6,
+				Eps1:       1e-8,
+				Eps2:       1e-8,
+			},
+			expected: []float64{1, 10, 1, 5, 4, 3},
+			tol:      1e-3,
+			settings: Settings{Iterations: 100, ObjectiveTol: 1e-16},
+		},
+		// Extended Rosenbrock problem.
+		LMTest{
+			prob: LMProblem{
+				Dim:        2,
+				Size:       2,
+				Func:       extendedRosenbrockFunc,
+				Jac:        rosenbrockNumJac.Jac,
+				InitParams: []float64{-20, 150},
+				Tau:        1e-6,
+				Eps1:       1e-8,
+				Eps2:       1e-8,
+			},
+			expected: []float64{1, 1},
+			tol:      1e-6,
+			settings: Settings{Iterations: 100, ObjectiveTol: 1e-16},
+		},
+	}
+
+	for _, testProb := range problems {
+		result, err := LM(testProb.prob, &testProb.settings)
+		if err != nil {
+			t.Errorf("unexepected error: %v", err)
+		}
+		if !floats.EqualApprox(result.X, testProb.expected, testProb.tol) {
+			t.Errorf("Optimal mismatch: got %v, want %v", result.X, testProb.expected)
+		}
+	}
+}
diff --git a/nlls/lmopt.go b/nlls/lmopt.go
new file mode 100644
index 0000000..6e0aac3
--- /dev/null
+++ b/nlls/lmopt.go
@@ -0,0 +1,223 @@
+// Copyright ©2019 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+Package nlls implements optimization routines for non-linear least squares problems
+using the Levenberg-Marquardt method. 
+
+Given function f:Rn -> Rm, where m is the number of non-linear functions and n parameters,
+the Levenberg-Marquardt method is used to seek a point X that minimizes F(x) = 0.5 * f.T * f. 
+
+The user supplies a non-linear function. The jacobian may also be supplied by the user or
+approximated by finite differences.
+*/
+package nlls
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/mat"
+	"gonum.org/v1/gonum/diff/fd"
+	"gonum.org/v1/gonum/optimize"
+)
+
+type Settings struct {
+	// Iterations represents the maximum number of iterations. Defaults to 100.
+	Iterations int
+
+	// ObjectiveTol represents the value for the obejective function after which
+	// the algorithm can stop. Defaults to 1e-16.
+	ObjectiveTol float64
+}
+
+func defaultSettings(set *Settings) {
+	set.Iterations = 100
+	set.ObjectiveTol = 1e-16
+}
+
+type Result struct {
+	X      []float64
+	Status optimize.Status
+}
+
+// NumJac is used if the user doesn't wish to provide a fucnction that evaluates 
+// the jacobian matrix. NumJac provides a method Jac that computes the jacobian matrix
+// by finite differences.
+type NumJac struct {
+	Func func(dst, param []float64)
+}
+
+func (nj *NumJac) Jac(dst *mat.Dense, param []float64) {
+	fd.Jacobian(dst, nj.Func, param, &fd.JacobianSettings{
+		Formula:    fd.Central,
+		Concurrent: true,
+	})
+}
+
+func maxDiagElem(m *mat.Dense) float64 {
+	r, c := m.Dims()
+	if r != c {
+		panic("lm: matrix is not square")
+	}
+	maxElem := m.At(0, 0)
+	for i := 1; i < r; i++ {
+		if m.At(i, i) > maxElem {
+			maxElem = m.At(i, i)
+		}
+	}
+	return maxElem
+}
+
+func addToDiag(m *mat.Dense, v float64) {
+	r, c := m.Dims()
+	if r != c {
+		panic("lm: matrix is not square")
+	}
+	for i := 0; i < r; i++ {
+		m.Set(i, i, m.At(i, i)+v)
+	}
+}
+
+func updateParams(dst []float64, params []float64, h *mat.VecDense) {
+	if len(params) != h.Len() {
+		panic("lm: lenghts don't match")
+	}
+	for i := 0; i < len(params); i++ {
+		dst[i] = params[i] - h.At(i, 0)
+	}
+}
+
+func calcRho(fParams []float64, fParamsNew []float64, h *mat.VecDense, grad *mat.VecDense, mu float64) float64 {
+	rho := floats.Dot(fParams, fParams) - floats.Dot(fParamsNew, fParamsNew)
+	tmpVec := mat.NewVecDense(h.Len(), nil)
+	tmpVec.AddScaledVec(grad, mu, h)
+	lDiff := mat.Dot(h, tmpVec)
+	rho /= lDiff
+	return rho
+}
+
+// LM is a function that solves non-linear least squares problems using the Levenberg-Marquardt
+// Method. 
+//
+// References:
+//  - Madsen, Kaj, Hans Bruun Nielsen, and Ole Tingleff. "Methods for non-linear least squares
+//    problems.", 2nd edition, 2004.
+//  - Lourakis, Manolis. "A Brief Description of the Levenberg-Marquardt Algorithm Implemened 
+//    by levmar", 2005.
+func LM(problem LMProblem, settings *Settings) (*Result, error) {
+	var set Settings
+	if settings != nil {
+		set = *settings
+	} else {
+		defaultSettings(&set)
+	}
+	dim := problem.Dim
+	if problem.Dim == 0 {
+		panic("lm: problem dimension is 0")
+	}
+	size := problem.Size
+	if problem.Size == 0 {
+		panic("lm: problem size is 0")
+	}
+	status := optimize.NotTerminated
+
+	dstFunc := make([]float64, size)
+	dstFuncNew := make([]float64, size)
+	dstJac := mat.NewDense(size, dim, nil)
+	dstA := mat.NewDense(dim, dim, nil)
+	dstGrad := mat.NewVecDense(dim, nil)
+	dstH := mat.NewVecDense(dim, nil)
+	nu := 2.0
+	var mu float64
+	found := false
+
+	// The inital guess is the zero vector by default.
+	parameters := make([]float64, dim)
+	parametersNew := make([]float64, dim)
+	if problem.InitParams != nil {
+		copy(parameters, problem.InitParams)
+	}
+
+	// Initial evaluation of A = J.T * J and g = J.T * f.
+	problem.Func(dstFunc, parameters)
+	problem.Jac(dstJac, parameters)
+	dstA.Mul(dstJac.T(), dstJac)
+	dstGrad.MulVec(dstJac.T(), mat.NewVecDense(size, dstFunc))
+
+	found = (mat.Norm(dstGrad, math.Inf(1)) <= problem.Eps1)
+	mu = problem.Tau * maxDiagElem(dstA)
+
+	for iter := 0; ; iter++ {
+		if iter == set.Iterations {
+			status = optimize.IterationLimit
+			break
+		}
+		if found {
+			status = optimize.StepConvergence
+			break
+		}
+
+		// Solve (A + mu * I) * h_lm = g.
+		addToDiag(dstA, mu)
+		err := dstH.SolveVec(dstA, dstGrad)
+		if err != nil {
+			panic("singular")
+		}
+
+		// Return A to its original state for the next steps. This is done in order not to copy A.
+		addToDiag(dstA, -mu)
+
+		if mat.Norm(dstH, 2) <= (floats.Norm(parameters, 2)+problem.Eps2)*problem.Eps2 {
+			found = true
+		} else {
+			updateParams(parametersNew, parameters, dstH)
+
+			// Calculate rho = (F(x) - F(x_new)) / (L(0) - L(h_lm)), where
+			// F = 0.5 * f.T * f, L = 0.5 * h_lm.T * (mu * h_lm - g).
+			problem.Func(dstFuncNew, parametersNew)
+			rho := calcRho(dstFunc, dstFuncNew, dstH, dstGrad, mu)
+
+			if rho > 0 { // step is acceptable
+				copy(parameters, parametersNew)
+				problem.Func(dstFunc, parameters)
+				problem.Jac(dstJac, parameters)
+				dstA.Mul(dstJac.T(), dstJac)
+				dstGrad.MulVec(dstJac.T(), mat.NewVecDense(size, dstFunc))
+				found = (mat.Norm(dstGrad, math.Inf(1)) <= problem.Eps1) ||
+					(0.5*floats.Dot(dstFunc, dstFunc) <= set.ObjectiveTol)
+				mu = mu * math.Max(1.0/3.0, 1-math.Pow(2*rho-1, 3))
+				nu = 2.0
+			} else {
+				mu *= nu
+				nu *= 2.0
+			}
+		}
+	}
+	return &Result{
+		X:      parameters,
+		Status: status,
+	}, nil
+}
+
+// LMProblem is used for running LM optimization. The objective function is
+// F = 0.5 * f.T * f, where f:Rn -> Rm and m >= n.
+type LMProblem struct {
+	// Dim is the dimension of the parameters of the problem (n).
+	Dim int
+	// Size specifies the number of nonlinear functions (m).
+	Size int
+	// Func computes the function value at params.
+	Func func(dst, param []float64)
+	// Jac computes the jacobian matrix of Func.
+	Jac func(dst *mat.Dense, param []float64)
+	// InitParams stores the users inital guess. Defaults to the zero vector when nil.
+	InitParams []float64
+	// Tau scales the initial damping parameter.
+	Tau float64
+	// Eps1 is a stopping criterion for the gradient of F.
+	Eps1 float64
+	// Eps2 is a stopping criterion for the step size.
+	Eps2 float64
+}