diff --git a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
index c4b19855a..007e13f9d 100644
--- a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
+++ b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
@@ -84,6 +84,12 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
 
     /** Current residuals. */
     protected double[] residuals;
+    
+    /** Weighted Jacobian */
+    protected double[][] wjacobian;
+    
+    /** Weighted residuals */
+    protected double[] wresiduals;
 
     /** Cost value (square root of the sum of the residuals). */
     protected double cost;
@@ -189,9 +195,10 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
         }
         for (int i = 0; i < rows; i++) {
             final double[] ji = jacobian[i];
-            final double factor = -Math.sqrt(residualsWeights[i]);
+            double wi = Math.sqrt(residualsWeights[i]);
             for (int j = 0; j < cols; ++j) {
-                ji[j] *= factor;
+                ji[j] *=  -1.0;
+                wjacobian[i][j] = ji[j]*wi;
             }
         }
     }
@@ -219,6 +226,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
         for (int i = 0; i < rows; i++) {
             final double residual = targetValues[i] - objective[i];
             residuals[i] = residual;
+            wresiduals[i]= residual*Math.sqrt(residualsWeights[i]);
             cost += residualsWeights[i] * residual * residual;
             index += cols;
         }
@@ -270,7 +278,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
             for (int j = i; j < cols; ++j) {
                 double sum = 0;
                 for (int k = 0; k < rows; ++k) {
-                    sum += jacobian[k][i] * jacobian[k][j];
+                    sum += wjacobian[k][i] * wjacobian[k][j];
                 }
                 jTj[i][j] = sum;
                 jTj[j][i] = sum;
@@ -342,6 +350,9 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
         cols      = point.length;
         jacobian  = new double[rows][cols];
 
+        wjacobian = new double[rows][cols];
+        wresiduals = new double[rows];
+        
         cost = Double.POSITIVE_INFINITY;
 
         return doOptimize();
diff --git a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
index b41456eff..286f82c7d 100644
--- a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
+++ b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
@@ -270,7 +270,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
         VectorialPointValuePair current = new VectorialPointValuePair(point, objective);
         while (true) {
             for (int i=0;i<rows;i++) {
-                qtf[i]=residuals[i];
+                qtf[i]=wresiduals[i];
             }
             incrementIterationsCounter();
 
@@ -285,7 +285,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             // so let jacobian contain the R matrix with its diagonal elements
             for (int k = 0; k < solvedCols; ++k) {
                 int pk = permutation[k];
-                jacobian[k][pk] = diagR[pk];
+                wjacobian[k][pk] = diagR[pk];
             }
 
             if (firstIteration) {
@@ -318,7 +318,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                     if (s != 0) {
                         double sum = 0;
                         for (int i = 0; i <= j; ++i) {
-                            sum += jacobian[i][pj] * qtf[i];
+                            sum += wjacobian[i][pj] * qtf[i];
                         }
                         maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
                     }
@@ -387,7 +387,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                     double dirJ = lmDir[pj];
                     work1[j] = 0;
                     for (int i = 0; i <= j; ++i) {
-                        work1[i] += jacobian[i][pj] * dirJ;
+                        work1[i] += wjacobian[i][pj] * dirJ;
                     }
                 }
                 double coeff1 = 0;
@@ -514,7 +514,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             int pk = permutation[k];
             double ypk = lmDir[pk] / diagR[pk];
             for (int i = 0; i < k; ++i) {
-                lmDir[permutation[i]] -= ypk * jacobian[i][pk];
+                lmDir[permutation[i]] -= ypk * wjacobian[i][pk];
             }
             lmDir[pk] = ypk;
         }
@@ -550,7 +550,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                 int pj = permutation[j];
                 double sum = 0;
                 for (int i = 0; i < j; ++i) {
-                    sum += jacobian[i][pj] * work1[permutation[i]];
+                    sum += wjacobian[i][pj] * work1[permutation[i]];
                 }
                 double s = (work1[pj] - sum) / diagR[pj];
                 work1[pj] = s;
@@ -565,7 +565,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             int pj = permutation[j];
             double sum = 0;
             for (int i = 0; i <= j; ++i) {
-                sum += jacobian[i][pj] * qy[i];
+                sum += wjacobian[i][pj] * qy[i];
             }
             sum /= diag[pj];
             sum2 += sum * sum;
@@ -625,7 +625,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                 work1[pj] /= work2[j];
                 double tmp = work1[pj];
                 for (int i = j + 1; i < solvedCols; ++i) {
-                    work1[permutation[i]] -= jacobian[i][pj] * tmp;
+                    work1[permutation[i]] -= wjacobian[i][pj] * tmp;
                 }
             }
             sum2 = 0;
@@ -676,7 +676,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
         for (int j = 0; j < solvedCols; ++j) {
             int pj = permutation[j];
             for (int i = j + 1; i < solvedCols; ++i) {
-                jacobian[i][pj] = jacobian[j][permutation[i]];
+                wjacobian[i][pj] = wjacobian[j][permutation[i]];
             }
             lmDir[j] = diagR[pj];
             work[j]  = qy[j];
@@ -707,7 +707,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
 
                     final double sin;
                     final double cos;
-                    double rkk = jacobian[k][pk];
+                    double rkk = wjacobian[k][pk];
                     if (Math.abs(rkk) < Math.abs(lmDiag[k])) {
                         final double cotan = rkk / lmDiag[k];
                         sin   = 1.0 / Math.sqrt(1.0 + cotan * cotan);
@@ -720,17 +720,17 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
 
                     // compute the modified diagonal element of R and
                     // the modified element of (Qty,0)
-                    jacobian[k][pk] = cos * rkk + sin * lmDiag[k];
+                    wjacobian[k][pk] = cos * rkk + sin * lmDiag[k];
                     final double temp = cos * work[k] + sin * qtbpj;
                     qtbpj = -sin * work[k] + cos * qtbpj;
                     work[k] = temp;
 
                     // accumulate the tranformation in the row of s
                     for (int i = k + 1; i < solvedCols; ++i) {
-                        double rik = jacobian[i][pk];
+                        double rik = wjacobian[i][pk];
                         final double temp2 = cos * rik + sin * lmDiag[i];
                         lmDiag[i] = -sin * rik + cos * lmDiag[i];
-                        jacobian[i][pk] = temp2;
+                        wjacobian[i][pk] = temp2;
                     }
 
                 }
@@ -738,8 +738,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
 
             // store the diagonal element of s and restore
             // the corresponding diagonal element of R
-            lmDiag[j] = jacobian[j][permutation[j]];
-            jacobian[j][permutation[j]] = lmDir[j];
+            lmDiag[j] = wjacobian[j][permutation[j]];
+            wjacobian[j][permutation[j]] = lmDir[j];
 
         }
 
@@ -759,7 +759,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                 int pj = permutation[j];
                 double sum = 0;
                 for (int i = j + 1; i < nSing; ++i) {
-                    sum += jacobian[i][pj] * work[i];
+                    sum += wjacobian[i][pj] * work[i];
                 }
                 work[j] = (work[j] - sum) / lmDiag[j];
             }
@@ -800,8 +800,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
         for (int k = 0; k < cols; ++k) {
             permutation[k] = k;
             double norm2 = 0;
-            for (int i = 0; i < jacobian.length; ++i) {
-                double akk = jacobian[i][k];
+            for (int i = 0; i < wjacobian.length; ++i) {
+                double akk = wjacobian[i][k];
                 norm2 += akk * akk;
             }
             jacNorm[k] = Math.sqrt(norm2);
@@ -815,8 +815,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             double ak2 = Double.NEGATIVE_INFINITY;
             for (int i = k; i < cols; ++i) {
                 double norm2 = 0;
-                for (int j = k; j < jacobian.length; ++j) {
-                    double aki = jacobian[j][permutation[i]];
+                for (int j = k; j < wjacobian.length; ++j) {
+                    double aki = wjacobian[j][permutation[i]];
                     norm2 += aki * aki;
                 }
                 if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
@@ -837,24 +837,24 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             permutation[k]          = pk;
 
             // choose alpha such that Hk.u = alpha ek
-            double akk   = jacobian[k][pk];
+            double akk   = wjacobian[k][pk];
             double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2);
             double betak = 1.0 / (ak2 - akk * alpha);
             beta[pk]     = betak;
 
             // transform the current column
             diagR[pk]        = alpha;
-            jacobian[k][pk] -= alpha;
+            wjacobian[k][pk] -= alpha;
 
             // transform the remaining columns
             for (int dk = cols - 1 - k; dk > 0; --dk) {
                 double gamma = 0;
-                for (int j = k; j < jacobian.length; ++j) {
-                    gamma += jacobian[j][pk] * jacobian[j][permutation[k + dk]];
+                for (int j = k; j < wjacobian.length; ++j) {
+                    gamma += wjacobian[j][pk] * wjacobian[j][permutation[k + dk]];
                 }
                 gamma *= betak;
-                for (int j = k; j < jacobian.length; ++j) {
-                    jacobian[j][permutation[k + dk]] -= gamma * jacobian[j][pk];
+                for (int j = k; j < wjacobian.length; ++j) {
+                    wjacobian[j][permutation[k + dk]] -= gamma * wjacobian[j][pk];
                 }
             }
 
@@ -874,11 +874,11 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             int pk = permutation[k];
             double gamma = 0;
             for (int i = k; i < rows; ++i) {
-                gamma += jacobian[i][pk] * y[i];
+                gamma += wjacobian[i][pk] * y[i];
             }
             gamma *= beta[pk];
             for (int i = k; i < rows; ++i) {
-                y[i] -= gamma * jacobian[i][pk];
+                y[i] -= gamma * wjacobian[i][pk];
             }
         }
     }
diff --git a/src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java b/src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
index ab70034bf..901e3150c 100644
--- a/src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
+++ b/src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
@@ -472,8 +472,8 @@ extends TestCase {
 
         VectorialPointValuePair optimum =
             optimizer.optimize(circle, target, weights, new double[] { 0, 0 });
-        assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1.0e-8);
-        assertEquals(0.2074999736353867,  optimum.getPointRef()[1], 1.0e-8);
+        assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1.0e-6);
+        assertEquals(0.2074999736353867,  optimum.getPointRef()[1], 1.0e-6);
         assertEquals(0.04268731682389561, optimizer.getRMS(),       1.0e-8);
 
     }