diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt
new file mode 100644
index 000000000..e6f575262
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticDifficultTest.kt
@@ -0,0 +1,90 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+
+fun main() {
+    val NData = 200
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1) - 104
+    }
+
+    val Nparams = 15
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat =  funcDifficultForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0 / Nparams * 1.0 - 0.085
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-2, 11.0, 9.0, 1.0)
+//    val opts = doubleArrayOf(3.0, 10000.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-3, 11.0, 9.0, 1.0)
+
+    val inputData = LMInput(::funcDifficultForLm,
+        p_init.as2D(),
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min.as2D(),
+        p_max.as2D(),
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        1)
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+    println("Y true and y received:")
+    var y_hat_after =  funcDifficultForLm(t_example, result.resultParameters, exampleNumber)
+    for (i in 0 until y_hat.shape.component1()) {
+        val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt
new file mode 100644
index 000000000..507943031
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticEasyTest.kt
@@ -0,0 +1,57 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcDifficultForLm
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcEasyForLm
+import space.kscience.kmath.tensors.LevenbergMarquardt.getStartDataForFuncEasy
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+
+fun main() {
+    val startedData = getStartDataForFuncEasy()
+    val inputData = LMInput(::funcEasyForLm,
+        DoubleTensorAlgebra.ones(ShapeND(intArrayOf(4, 1))).as2D(),
+        startedData.t,
+        startedData.y_dat,
+        startedData.weight,
+        startedData.dp,
+        startedData.p_min,
+        startedData.p_max,
+        startedData.opts[1].toInt(),
+        doubleArrayOf(startedData.opts[2], startedData.opts[3], startedData.opts[4], startedData.opts[5]),
+        doubleArrayOf(startedData.opts[6], startedData.opts[7], startedData.opts[8]),
+        startedData.opts[9].toInt(),
+        10,
+        startedData.example_number)
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+    println("Y true and y received:")
+    var y_hat_after =  funcDifficultForLm(startedData.t, result.resultParameters, startedData.example_number)
+    for (i in 0 until startedData.y_dat.shape.component1()) {
+        val x = (startedData.y_dat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt
new file mode 100644
index 000000000..0659db103
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StaticLm/staticMiddleTest.kt
@@ -0,0 +1,88 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt.StaticLm
+
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.LevenbergMarquardt.funcMiddleForLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.math.roundToInt
+fun main() {
+    val NData = 100
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1)
+    }
+
+    val Nparams = 20
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat =  funcMiddleForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0)
+
+    val inputData = LMInput(::funcMiddleForLm,
+        p_init.as2D(),
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min.as2D(),
+        p_max.as2D(),
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        1)
+
+    val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+    println("Parameters:")
+    for (i in 0 until result.resultParameters.shape.component1()) {
+        val x = (result.resultParameters[i, 0] * 10000).roundToInt() / 10000.0
+        print("$x ")
+    }
+    println()
+
+
+    var y_hat_after =  funcMiddleForLm(t_example, result.resultParameters, exampleNumber)
+    for (i in 0 until y_hat.shape.component1()) {
+        val x = (y_hat[i, 0] * 10000).roundToInt() / 10000.0
+        val y = (y_hat_after[i, 0] * 10000).roundToInt() / 10000.0
+        println("$x $y")
+    }
+
+    println("Сhi_sq:")
+    println(result.resultChiSq)
+    println("Number of iterations:")
+    println(result.iterations)
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
new file mode 100644
index 000000000..b2818ef2a
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
@@ -0,0 +1,68 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt.StreamingLm
+
+import kotlinx.coroutines.delay
+import kotlinx.coroutines.flow.*
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.tensors.LevenbergMarquardt.StartDataLm
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.LMInput
+import space.kscience.kmath.tensors.core.levenbergMarquardt
+import kotlin.random.Random
+import kotlin.reflect.KFunction3
+
+fun streamLm(lm_func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
+             startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int): Flow<MutableStructure2D<Double>> = flow{
+
+    var example_number = startData.example_number
+    var p_init = startData.p_init
+    var t = startData.t
+    var y_dat = startData.y_dat
+    val weight = startData.weight
+    val dp = startData.dp
+    val p_min = startData.p_min
+    val p_max = startData.p_max
+    val opts = startData.opts
+
+    var steps = numberOfLaunches
+    val isEndless = (steps <= 0)
+
+    val inputData = LMInput(lm_func,
+        p_init,
+        t,
+        y_dat,
+        weight,
+        dp,
+        p_min,
+        p_max,
+        opts[1].toInt(),
+        doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+        doubleArrayOf(opts[6], opts[7], opts[8]),
+        opts[9].toInt(),
+        10,
+        example_number)
+
+    while (isEndless || steps > 0) {
+        val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+        emit(result.resultParameters)
+        delay(launchFrequencyInMs)
+        inputData.realValues = generateNewYDat(y_dat, 0.1)
+        inputData.startParameters = result.resultParameters
+        if (!isEndless) steps -= 1
+    }
+}
+
+fun generateNewYDat(y_dat: MutableStructure2D<Double>, delta: Double): MutableStructure2D<Double>{
+    val n = y_dat.shape.component1()
+    val y_dat_new = zeros(ShapeND(intArrayOf(n, 1))).as2D()
+    for (i in 0 until n) {
+        val randomEps = Random.nextDouble(delta + delta) - delta
+        y_dat_new[i, 0] = y_dat[i, 0] + randomEps
+    }
+    return y_dat_new
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt
new file mode 100644
index 000000000..c9dd5029e
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamingLmTest.kt
@@ -0,0 +1,32 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt.StreamingLm
+
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.tensors.LevenbergMarquardt.*
+import kotlin.math.roundToInt
+
+suspend fun main(){
+    val startData = getStartDataForFuncDifficult()
+    // Создание потока:
+    val lmFlow = streamLm(::funcDifficultForLm, startData, 0,  100)
+    var initialTime = System.currentTimeMillis()
+    var lastTime: Long
+    val launches = mutableListOf<Long>()
+    // Запуск потока
+    lmFlow.collect { parameters ->
+        lastTime = System.currentTimeMillis()
+        launches.add(lastTime - initialTime)
+        initialTime = lastTime
+        for (i in 0 until parameters.shape.component1()) {
+            val x = (parameters[i, 0] * 10000).roundToInt() / 10000.0
+            print("$x ")
+            if (i == parameters.shape.component1() - 1) println()
+        }
+    }
+
+    println("Average without first is: ${launches.subList(1, launches.size - 1).average()}")
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
new file mode 100644
index 000000000..7ccb37ed0
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
@@ -0,0 +1,222 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.LevenbergMarquardt
+
+import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.pow
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times
+import space.kscience.kmath.tensors.core.asDoubleTensor
+
+public data class StartDataLm (
+    var lm_matx_y_dat: MutableStructure2D<Double>,
+    var example_number: Int,
+    var p_init: MutableStructure2D<Double>,
+    var t: MutableStructure2D<Double>,
+    var y_dat: MutableStructure2D<Double>,
+    var weight: Double,
+    var dp: MutableStructure2D<Double>,
+    var p_min: MutableStructure2D<Double>,
+    var p_max: MutableStructure2D<Double>,
+    var consts: MutableStructure2D<Double>,
+    var opts: DoubleArray
+)
+
+fun funcEasyForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))
+
+    if (exampleNumber == 1) {
+        y_hat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))).times(p[0, 0]) + t.times(p[2, 0]).times(
+            DoubleTensorAlgebra.exp((t.times(-1.0 / p[3, 0])))
+        )
+    }
+    else if (exampleNumber == 2) {
+        val mt = t.max()
+        y_hat = (t.times(1.0 / mt)).times(p[0, 0]) +
+                (t.times(1.0 / mt)).pow(2).times(p[1, 0]) +
+                (t.times(1.0 / mt)).pow(3).times(p[2, 0]) +
+                (t.times(1.0 / mt)).pow(4).times(p[3, 0])
+    }
+    else if (exampleNumber == 3) {
+        y_hat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0])))
+            .times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0])
+    }
+
+    return y_hat.as2D()
+}
+
+fun funcMiddleForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1)))
+
+    val mt = t.max()
+    for(i in 0 until p.shape.component1()){
+        y_hat += (t.times(1.0 / mt)).times(p[i, 0])
+    }
+
+    for(i in 0 until 5){
+        y_hat = funcEasyForLm(y_hat.as2D(), p, exampleNumber).asDoubleTensor()
+    }
+
+    return y_hat.as2D()
+}
+
+fun funcDifficultForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+    val m = t.shape.component1()
+    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1)))
+
+    val mt = t.max()
+    for(i in 0 until p.shape.component1()){
+        y_hat = y_hat.plus( (t.times(1.0 / mt)).times(p[i, 0]) )
+    }
+
+    for(i in 0 until 4){
+        y_hat = funcEasyForLm((y_hat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor()
+    }
+
+    return y_hat.as2D()
+}
+
+
+fun getStartDataForFuncDifficult(): StartDataLm  {
+    val NData = 200
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1) - 104
+    }
+
+    val Nparams = 15
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat =  funcDifficultForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 0.9)
+    }
+
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0 / Nparams * 1.0 - 0.085
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0)
+
+    return StartDataLm(y_dat, 1, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts)
+}
+
+fun getStartDataForFuncMiddle(): StartDataLm  {
+    val NData = 100
+    var t_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+    for (i in 0 until NData) {
+        t_example[i, 0] = t_example[i, 0] * (i + 1)
+    }
+
+    val Nparams = 20
+    var p_example = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_example[i, 0] = p_example[i, 0] + i - 25
+    }
+
+    val exampleNumber = 1
+
+    var y_hat =  funcMiddleForLm(t_example, p_example, exampleNumber)
+
+    var p_init = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+    for (i in 0 until Nparams) {
+        p_init[i, 0] = (p_example[i, 0] + 10.0)
+    }
+    var t = t_example
+    val y_dat = y_hat
+    val weight = 1.0
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+    var p_min = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / -50.0)
+    val p_max = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+    p_min = p_min.div(1.0 / 50.0)
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+    val opts = doubleArrayOf(3.0, 10000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0)
+
+    var example_number = 1
+
+    return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min.as2D(), p_max.as2D(), consts, opts)
+}
+
+fun getStartDataForFuncEasy(): StartDataLm {
+    val lm_matx_y_dat = doubleArrayOf(
+        19.6594, 18.6096, 17.6792, 17.2747, 16.3065, 17.1458, 16.0467, 16.7023, 15.7809, 15.9807,
+        14.7620, 15.1128, 16.0973, 15.1934, 15.8636, 15.4763, 15.6860, 15.1895, 15.3495, 16.6054,
+        16.2247, 15.9854, 16.1421, 17.0960, 16.7769, 17.1997, 17.2767, 17.5882, 17.5378, 16.7894,
+        17.7648, 18.2512, 18.1581, 16.7037, 17.8475, 17.9081, 18.3067, 17.9632, 18.2817, 19.1427,
+        18.8130, 18.5658, 18.0056, 18.4607, 18.5918, 18.2544, 18.3731, 18.7511, 19.3181, 17.3066,
+        17.9632, 19.0513, 18.7528, 18.2928, 18.5967, 17.8567, 17.7859, 18.4016, 18.9423, 18.4959,
+        17.8000, 18.4251, 17.7829, 17.4645, 17.5221, 17.3517, 17.4637, 17.7563, 16.8471, 17.4558,
+        17.7447, 17.1487, 17.3183, 16.8312, 17.7551, 17.0942, 15.6093, 16.4163, 15.3755, 16.6725,
+        16.2332, 16.2316, 16.2236, 16.5361, 15.3721, 15.3347, 15.5815, 15.6319, 14.4538, 14.6044,
+        14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157
+    )
+
+    var example_number = 1
+    val p_init = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0)
+    ).as2D()
+
+    var t = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(100, 1))).as2D()
+    for (i in 0 until 100) {
+        t[i, 0] = t[i, 0] * (i + 1)
+    }
+
+    val y_dat = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(100, 1)), lm_matx_y_dat
+    ).as2D()
+
+    val weight = 4.0
+
+    val dp = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+    ).as2D()
+
+    val p_min = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(-50.0, -20.0, -2.0, -100.0)
+    ).as2D()
+
+    val p_max = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0)
+    ).as2D()
+
+    val consts = BroadcastDoubleTensorAlgebra.fromArray(
+        ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.0)
+    ).as2D()
+
+    val opts = doubleArrayOf(3.0, 100.0, 1e-3, 1e-3, 1e-1, 1e-1, 1e-2, 11.0, 9.0, 1.0)
+
+    return StartDataLm(y_dat, example_number, p_init, t, y_dat, weight, dp, p_min, p_max, consts, opts)
+}
\ No newline at end of file
diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
index c56583fa8..8ad7f7293 100644
--- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
@@ -19,9 +19,9 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
 import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
-import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.FloatField
diff --git a/kmath-symja/build.gradle.kts b/kmath-symja/build.gradle.kts
index 8741de2ae..a996f3bec 100644
--- a/kmath-symja/build.gradle.kts
+++ b/kmath-symja/build.gradle.kts
@@ -10,7 +10,7 @@ plugins {
 description = "Symja integration module"
 
 dependencies {
-    api("org.matheclipse:matheclipse-core:2.0.0-SNAPSHOT") {
+    api("org.matheclipse:matheclipse-core:2.0.0") {
         // Incorrect transitive dependencies
         exclude("org.apfloat", "apfloat")
         exclude("org.hipparchus", "hipparchus-clustering")
diff --git a/kmath-tensors/build.gradle.kts b/kmath-tensors/build.gradle.kts
index 79c39bae7..2497314f0 100644
--- a/kmath-tensors/build.gradle.kts
+++ b/kmath-tensors/build.gradle.kts
@@ -4,7 +4,13 @@ plugins {
 
 kscience{
     jvm()
-    js()
+    js {
+        browser {
+            testTask {
+                useMocha().timeout = "0"
+            }
+        }
+    }
     native()
 
     dependencies {
diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt
index faff2eb80..f2c7f1821 100644
--- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/LinearOpsTensorAlgebra.kt
@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.tensors.api
 
+import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.Field
 
@@ -103,4 +104,11 @@ public interface LinearOpsTensorAlgebra<T, A : Field<T>> : TensorPartialDivision
      */
     public fun symEig(structureND: StructureND<T>): Pair<StructureND<T>, StructureND<T>>
 
+    /** Returns the solution to the equation Ax = B for the square matrix A as `input1` and
+     * for the square matrix B as `input2`.
+     *
+     * @receiver the `input1` and the `input2`.
+     * @return the square matrix x which is the solution of the equation.
+     */
+    public fun solve(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double>
 }
diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
index 70a3ef7e2..c8cf56888 100644
--- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
@@ -17,10 +17,7 @@ import space.kscience.kmath.tensors.api.AnalyticTensorAlgebra
 import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra
 import space.kscience.kmath.tensors.api.Tensor
 import space.kscience.kmath.tensors.core.internal.*
-import kotlin.math.abs
-import kotlin.math.ceil
-import kotlin.math.floor
-import kotlin.math.sqrt
+import kotlin.math.*
 
 /**
  * Implementation of basic operations over double tensors and basic algebra operations on them.
@@ -706,14 +703,18 @@ public open class DoubleTensorAlgebra :
     override fun svd(
         structureND: StructureND<Double>,
     ): Triple<StructureND<Double>, StructureND<Double>, StructureND<Double>> =
-        svd(structureND = structureND, epsilon = 1e-10)
+        svdGolubKahan(structureND = structureND, epsilon = 1e-10)
 
     override fun symEig(structureND: StructureND<Double>): Pair<DoubleTensor, DoubleTensor> =
         symEigJacobi(structureND = structureND, maxIteration = 50, epsilon = 1e-15)
 
+    override fun solve(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
+        val aSvd = DoubleTensorAlgebra.svd(a)
+        val s = BroadcastDoubleTensorAlgebra.diagonalEmbedding(aSvd.second.map {1.0 / it})
+        val aInverse = aSvd.third.dot(s).dot(aSvd.first.transposed())
+        return aInverse.dot(b).as2D()
+    }
 }
 
 public val Double.Companion.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra
-public val DoubleField.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra
-
-
+public val DoubleField.tensorAlgebra: DoubleTensorAlgebra get() = DoubleTensorAlgebra
\ No newline at end of file
diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt
new file mode 100644
index 000000000..3cb485d7d
--- /dev/null
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt
@@ -0,0 +1,589 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.core
+
+import space.kscience.kmath.linear.transpose
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.dot
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.minus
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.times
+import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.transposed
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus
+import kotlin.math.max
+import kotlin.math.min
+import kotlin.math.pow
+import kotlin.reflect.KFunction3
+
+/**
+ * Type of convergence achieved as a result of executing the Levenberg-Marquardt algorithm.
+ *
+ * InGradient: gradient convergence achieved
+ *            (max(J^T W dy) < epsilon1,
+ *            where J - Jacobi matrix (dy^/dp) for the current approximation y^,
+ *            W - weight matrix from input, dy = (y - y^(p))).
+ * InParameters: convergence in parameters achieved
+ *            (max(h_i / p_i) < epsilon2,
+ *            where h_i - offset for parameter p_i on the current iteration).
+ * InReducedChiSquare: chi-squared convergence achieved
+ *            (chi squared value divided by (m - n + 1) < epsilon2,
+ *            where n - number of parameters, m - amount of points).
+ * NoConvergence: the maximum number of iterations has been reached without reaching any convergence.
+ */
+public enum class TypeOfConvergence {
+    InGradient,
+    InParameters,
+    InReducedChiSquare,
+    NoConvergence
+}
+
+/**
+ * The data obtained as a result of the execution of the Levenberg-Marquardt algorithm.
+ *
+ *  iterations: number of completed iterations.
+ *  funcCalls: the number of evaluations of the input function during execution.
+ *  resultChiSq: chi squared value on final parameters.
+ *  resultLambda: final lambda parameter used to calculate the offset.
+ *  resultParameters: final parameters.
+ *  typeOfConvergence: type of convergence.
+ */
+public data class LMResultInfo (
+    var iterations:Int,
+    var funcCalls: Int,
+    var resultChiSq: Double,
+    var resultLambda: Double,
+    var resultParameters: MutableStructure2D<Double>,
+    var typeOfConvergence: TypeOfConvergence,
+)
+
+/**
+ * Input data for the Levenberg-Marquardt function.
+ *
+ *  func: function of n independent variables x, m parameters an example number,
+ *        rotating a vector of n values y, in which each of the y_i is calculated at its x_i with the given parameters.
+ *  startParameters: starting parameters.
+ *  independentVariables: independent variables, for each of which the real value is known.
+ *  realValues: real values obtained with given independent variables but unknown parameters.
+ *  weight: measurement error for realValues (denominator in each term of sum of weighted squared errors).
+ *  pDelta: delta when calculating the derivative with respect to parameters.
+ *  minParameters: the lower bound of parameter values.
+ *  maxParameters: upper limit of parameter values.
+ *  maxIterations: maximum allowable number of iterations.
+ *  epsilons: epsilon1 - convergence tolerance for gradient,
+ *            epsilon2 - convergence tolerance for parameters,
+ *            epsilon3 - convergence tolerance for reduced chi-square,
+ *            epsilon4 - determines acceptance of a step.
+ *  lambdas: lambda0 - starting lambda value for parameter offset count,
+ *           lambdaUp - factor for increasing lambda,
+ *           lambdaDown - factor for decreasing lambda.
+ *  updateType: 1: Levenberg-Marquardt lambda update,
+ *              2: Quadratic update,
+ *              3: Nielsen's lambda update equations.
+ *  nargin: a value that determines which options to use by default
+ *         (<5 - use weight by default, <6 - use pDelta by default, <7 - use minParameters by default,
+ *          <8 - use maxParameters by default, <9 - use updateType by default).
+ *  exampleNumber: a parameter for a function with which you can choose its behavior.
+ */
+public data class LMInput (
+    var func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
+    var startParameters: MutableStructure2D<Double>,
+    var independentVariables: MutableStructure2D<Double>,
+    var realValues: MutableStructure2D<Double>,
+    var weight: Double,
+    var pDelta: MutableStructure2D<Double>,
+    var minParameters: MutableStructure2D<Double>,
+    var maxParameters: MutableStructure2D<Double>,
+    var maxIterations: Int,
+    var epsilons: DoubleArray,
+    var lambdas: DoubleArray,
+    var updateType: Int,
+    var nargin: Int,
+    var exampleNumber: Int
+)
+
+
+/**
+ * Levenberg-Marquardt optimization.
+ *
+ * An optimization method that iteratively searches for the optimal function parameters
+ * that best describe the dataset. The 'input' is the function being optimized, a set of real data
+ * (calculated with independent variables, but with an unknown set of parameters), a set of
+ * independent variables, and variables for adjusting the algorithm, described in the documentation for the LMInput class.
+ * The function returns number of completed iterations, the number of evaluations of the input function during execution,
+ * chi squared value on final parameters, final lambda parameter used to calculate the offset, final parameters
+ * and type of convergence in the 'output'.
+ *
+ * @receiver the `input`.
+ * @return the 'output'.
+ */
+public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultInfo {
+    val resultInfo = LMResultInfo(0, 0, 0.0,
+        0.0, inputData.startParameters, TypeOfConvergence.NoConvergence)
+
+    val eps = 2.2204e-16
+
+    val settings = LMSettings(0, 0, inputData.exampleNumber)
+    settings.funcCalls = 0 // running count of function evaluations
+
+    var p = inputData.startParameters
+    val t = inputData.independentVariables
+
+    val Npar   = length(p) // number of parameters
+    val Npnt   = length(inputData.realValues) // number of data points
+    var pOld = zeros(ShapeND(intArrayOf(Npar, 1))).as2D() // previous set of parameters
+    var yOld = zeros(ShapeND(intArrayOf(Npnt, 1))).as2D() // previous model, y_old = y_hat(t;p_old)
+    var X2 = 1e-3 / eps  // a really big initial Chi-sq value
+    var X2Old = 1e-3 / eps // a really big initial Chi-sq value
+    var J = zeros(ShapeND(intArrayOf(Npnt, Npar))).as2D() // Jacobian matrix
+    val DoF = Npnt - Npar // statistical degrees of freedom
+
+    var weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(inputData.weight)).as2D()
+    if (inputData.nargin <  5) {
+        weight = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])).as2D()
+    }
+
+    var dp = inputData.pDelta
+    if (inputData.nargin < 6) {
+        dp = fromArray(ShapeND(intArrayOf(1, 1)), doubleArrayOf(0.001)).as2D()
+    }
+
+    var minParameters = inputData.minParameters
+    if (inputData.nargin < 7) {
+        minParameters = p
+        minParameters.abs()
+        minParameters = minParameters.div(-100.0).as2D()
+    }
+
+    var maxParameters = inputData.maxParameters
+    if (inputData.nargin < 8) {
+        maxParameters = p
+        maxParameters.abs()
+        maxParameters = maxParameters.div(100.0).as2D()
+    }
+
+    var maxIterations = inputData.maxIterations
+    var epsilon1 = inputData.epsilons[0]
+    var epsilon2 = inputData.epsilons[1]
+    var epsilon3 = inputData.epsilons[2]
+    var epsilon4 = inputData.epsilons[3]
+    var lambda0  = inputData.lambdas[0]
+    var lambdaUpFac = inputData.lambdas[1]
+    var lambdaDnFac = inputData.lambdas[2]
+    var updateType = inputData.updateType
+
+    if (inputData.nargin < 9) {
+        maxIterations = 10 * Npar
+        epsilon1 = 1e-3
+        epsilon2 = 1e-3
+        epsilon3 = 1e-1
+        epsilon4 = 1e-1
+        lambda0 = 1e-2
+        lambdaUpFac = 11.0
+        lambdaDnFac = 9.0
+        updateType = 1
+    }
+
+    minParameters = makeColumn(minParameters)
+    maxParameters = makeColumn(maxParameters)
+
+    if (length(makeColumn(dp)) == 1) {
+        dp = ones(ShapeND(intArrayOf(Npar, 1))).div(1 / dp[0, 0]).as2D()
+    }
+
+    var stop = false // termination flag
+
+    if (weight.shape.component1() == 1 || variance(weight) == 0.0) { // identical weights vector
+        weight = ones(ShapeND(intArrayOf(Npnt, 1))).div(1 / kotlin.math.abs(weight[0, 0])).as2D()
+    }
+    else {
+        weight = makeColumn(weight)
+        weight.abs()
+    }
+
+    // initialize Jacobian with finite difference calculation
+    var lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, 1, J, p, inputData.realValues, weight, dp, settings)
+    var JtWJ = lmMatxAns[0]
+    var JtWdy = lmMatxAns[1]
+    X2 = lmMatxAns[2][0, 0]
+    var yHat = lmMatxAns[3]
+    J = lmMatxAns[4]
+
+    if ( abs(JtWdy).max() < epsilon1 ) {
+        stop = true
+    }
+
+    var lambda = 1.0
+    var nu = 1
+
+    if (updateType == 1) {
+        lambda  = lambda0 // Marquardt: init'l lambda
+    }
+    else {
+        lambda  = lambda0 * (makeColumnFromDiagonal(JtWJ)).max()
+        nu = 2
+    }
+
+    X2Old = X2 // previous value of X2
+
+    var h: DoubleTensor
+
+    while (!stop && settings.iteration <= maxIterations) {
+        settings.iteration += 1
+
+        // incremental change in parameters
+        h = if (updateType == 1) { // Marquardt
+            val solve = solve(JtWJ.plus(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda)).as2D(), JtWdy)
+            solve.asDoubleTensor()
+        } else { // Quadratic and Nielsen
+            val solve = solve(JtWJ.plus(lmEye(Npar).div(1 / lambda)).as2D(), JtWdy)
+            solve.asDoubleTensor()
+        }
+
+        var pTry = (p + h).as2D()  // update the [idx] elements
+        pTry = smallestElementComparison(largestElementComparison(minParameters, pTry.as2D()), maxParameters) // apply constraints
+
+        var deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try
+
+        for (i in 0 until deltaY.shape.component1()) {  // floating point error; break
+            for (j in 0 until deltaY.shape.component2()) {
+                if (deltaY[i, j] == Double.POSITIVE_INFINITY || deltaY[i, j] == Double.NEGATIVE_INFINITY) {
+                    stop = true
+                    break
+                }
+            }
+        }
+
+        settings.funcCalls += 1
+
+        val tmp = deltaY.times(weight)
+        var X2Try = deltaY.as2D().transpose().dot(tmp) // Chi-squared error criteria
+
+        val alpha = 1.0
+        if (updateType == 2) { // Quadratic
+            // One step of quadratic line update in the h direction for minimum X2
+            val alpha = JtWdy.transpose().dot(h) / ((X2Try.minus(X2)).div(2.0).plus(2 * JtWdy.transpose().dot(h)))
+            h = h.dot(alpha)
+            pTry = p.plus(h).as2D() // update only [idx] elements
+            pTry = smallestElementComparison(largestElementComparison(minParameters, pTry), maxParameters) // apply constraints
+
+            deltaY = inputData.realValues.minus(evaluateFunction(inputData.func, t, pTry, inputData.exampleNumber)) // residual error using p_try
+            settings.funcCalls += 1
+
+            X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria
+        }
+
+        val rho = when (updateType) { // Nielsen
+            1 -> {
+                val tmp = h.transposed()
+                    .dot(makeMatrixWithDiagonal(makeColumnFromDiagonal(JtWJ)).div(1 / lambda).dot(h).plus(JtWdy))
+                X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0]
+            }
+            else -> {
+                val tmp = h.transposed().dot(h.div(1 / lambda).plus(JtWdy))
+                X2.minus(X2Try).as2D()[0, 0] / abs(tmp.as2D()).as2D()[0, 0]
+            }
+        }
+
+        if (rho > epsilon4) { // it IS significantly better
+            val dX2 = X2.minus(X2Old)
+            X2Old = X2
+            pOld = p.copyToTensor().as2D()
+            yOld = yHat.copyToTensor().as2D()
+            p = makeColumn(pTry) // accept p_try
+
+            lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, dX2.toInt(), J, p, inputData.realValues, weight, dp, settings)
+            // decrease lambda ==> Gauss-Newton method
+            JtWJ = lmMatxAns[0]
+            JtWdy = lmMatxAns[1]
+            X2 = lmMatxAns[2][0, 0]
+            yHat = lmMatxAns[3]
+            J = lmMatxAns[4]
+
+            lambda = when (updateType) {
+                1 -> { // Levenberg
+                    max(lambda / lambdaDnFac, 1e-7);
+                }
+
+                2 -> { // Quadratic
+                    max(lambda / (1 + alpha), 1e-7);
+                }
+
+                else -> { // Nielsen
+                    nu = 2
+                    lambda * max(1.0 / 3, 1 - (2 * rho - 1).pow(3))
+                }
+            }
+        } else { // it IS NOT better
+            X2 = X2Old // do not accept p_try
+            if (settings.iteration % (2 * Npar) == 0) { // rank-1 update of Jacobian
+                lmMatxAns = lmMatx(inputData.func, t, pOld, yOld, -1, J, p, inputData.realValues, weight, dp, settings)
+                JtWJ = lmMatxAns[0]
+                JtWdy = lmMatxAns[1]
+                yHat = lmMatxAns[3]
+                J = lmMatxAns[4]
+            }
+
+            // increase lambda  ==> gradient descent method
+            lambda = when (updateType) {
+                1 -> { // Levenberg
+                    min(lambda * lambdaUpFac, 1e7)
+                }
+
+                2 -> { // Quadratic
+                    lambda + kotlin.math.abs(((X2Try.as2D()[0, 0] - X2) / 2) / alpha)
+                }
+
+                else -> { // Nielsen
+                    nu *= 2
+                    lambda * (nu / 2)
+                }
+            }
+        }
+
+        val chiSq = X2 / DoF
+        resultInfo.iterations = settings.iteration
+        resultInfo.funcCalls = settings.funcCalls
+        resultInfo.resultChiSq = chiSq
+        resultInfo.resultLambda = lambda
+        resultInfo.resultParameters = p
+
+
+        if (abs(JtWdy).max() < epsilon1 && settings.iteration > 2) {
+            resultInfo.typeOfConvergence = TypeOfConvergence.InGradient
+            stop = true
+        }
+        if ((abs(h.as2D()).div(abs(p) + 1e-12)).max() < epsilon2 && settings.iteration > 2) {
+            resultInfo.typeOfConvergence = TypeOfConvergence.InParameters
+            stop = true
+        }
+        if (X2 / DoF < epsilon3 && settings.iteration > 2) {
+            resultInfo.typeOfConvergence = TypeOfConvergence.InReducedChiSquare
+            stop = true
+        }
+        if (settings.iteration == maxIterations) {
+            resultInfo.typeOfConvergence = TypeOfConvergence.NoConvergence
+            stop = true
+        }
+    }
+    return resultInfo
+}
+
+private data class LMSettings (
+    var iteration:Int,
+    var funcCalls: Int,
+    var exampleNumber:Int
+)
+
+/* matrix -> column of all elements */
+private fun makeColumn(tensor: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val shape = intArrayOf(tensor.shape.component1() * tensor.shape.component2(), 1)
+    val buffer = DoubleArray(tensor.shape.component1() * tensor.shape.component2())
+    for (i in 0 until tensor.shape.component1()) {
+        for (j in 0 until tensor.shape.component2()) {
+            buffer[i * tensor.shape.component2() + j] = tensor[i, j]
+        }
+    }
+    return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(shape), buffer).as2D()
+}
+
+/* column length */
+private fun length(column: MutableStructure2D<Double>) : Int {
+    return column.shape.component1()
+}
+
+private fun MutableStructure2D<Double>.abs() {
+    for (i in 0 until this.shape.component1()) {
+        for (j in 0 until this.shape.component2()) {
+            this[i, j] = kotlin.math.abs(this[i, j])
+        }
+    }
+}
+
+private fun abs(input: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val tensor = BroadcastDoubleTensorAlgebra.ones(
+        ShapeND(
+            intArrayOf(
+                input.shape.component1(),
+                input.shape.component2()
+            )
+        )
+    ).as2D()
+    for (i in 0 until tensor.shape.component1()) {
+        for (j in 0 until tensor.shape.component2()) {
+            tensor[i, j] = kotlin.math.abs(input[i, j])
+        }
+    }
+    return tensor
+}
+
+private fun makeColumnFromDiagonal(input: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val tensor = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(input.shape.component1(), 1))).as2D()
+    for (i in 0 until tensor.shape.component1()) {
+        tensor[i, 0] = input[i, i]
+    }
+    return tensor
+}
+
+private fun makeMatrixWithDiagonal(column: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val size = column.shape.component1()
+    val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(size, size))).as2D()
+    for (i in 0 until size) {
+        tensor[i, i] = column[i, 0]
+    }
+    return tensor
+}
+
+private fun lmEye(size: Int): MutableStructure2D<Double> {
+    val column = BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(size, 1))).as2D()
+    return makeMatrixWithDiagonal(column)
+}
+
+private fun largestElementComparison(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val aSizeX = a.shape.component1()
+    val aSizeY = a.shape.component2()
+    val bSizeX = b.shape.component1()
+    val bSizeY = b.shape.component2()
+    val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D()
+    for (i in 0 until tensor.shape.component1()) {
+        for (j in 0 until tensor.shape.component2()) {
+            if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) {
+                tensor[i, j] = max(a[i, j], b[i, j])
+            }
+            else if (i < aSizeX && j < aSizeY) {
+                tensor[i, j] = a[i, j]
+            }
+            else {
+                tensor[i, j] = b[i, j]
+            }
+        }
+    }
+    return tensor
+}
+
+private fun smallestElementComparison(a: MutableStructure2D<Double>, b: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    val aSizeX = a.shape.component1()
+    val aSizeY = a.shape.component2()
+    val bSizeX = b.shape.component1()
+    val bSizeY = b.shape.component2()
+    val tensor = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(max(aSizeX, bSizeX), max(aSizeY, bSizeY)))).as2D()
+    for (i in 0 until tensor.shape.component1()) {
+        for (j in 0 until tensor.shape.component2()) {
+            if (i < aSizeX && i < bSizeX && j < aSizeY && j < bSizeY) {
+                tensor[i, j] = min(a[i, j], b[i, j])
+            }
+            else if (i < aSizeX && j < aSizeY) {
+                tensor[i, j] = a[i, j]
+            }
+            else {
+                tensor[i, j] = b[i, j]
+            }
+        }
+    }
+    return tensor
+}
+
+private fun getZeroIndices(column: MutableStructure2D<Double>, epsilon: Double = 0.000001): MutableStructure2D<Double>? {
+    var idx = emptyArray<Double>()
+    for (i in 0 until column.shape.component1()) {
+        if (kotlin.math.abs(column[i, 0]) > epsilon) {
+            idx += (i + 1.0)
+        }
+    }
+    if (idx.isNotEmpty()) {
+        return BroadcastDoubleTensorAlgebra.fromArray(ShapeND(intArrayOf(idx.size, 1)), idx.toDoubleArray()).as2D()
+    }
+    return null
+}
+
+private fun evaluateFunction(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) ->  MutableStructure2D<Double>,
+                             t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int)
+        : MutableStructure2D<Double>
+{
+    return func(t, p, exampleNumber)
+}
+
+private fun lmMatx(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> MutableStructure2D<Double>,
+                   t: MutableStructure2D<Double>, pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>,
+                   dX2: Int, JInput: MutableStructure2D<Double>, p: MutableStructure2D<Double>,
+                   yDat: MutableStructure2D<Double>, weight: MutableStructure2D<Double>, dp:MutableStructure2D<Double>, settings:LMSettings) : Array<MutableStructure2D<Double>>
+{
+    // default: dp = 0.001
+    val Npar = length(p) // number of parameters
+
+    val yHat = evaluateFunction(func, t, p, settings.exampleNumber) // evaluate model using parameters 'p'
+    settings.funcCalls += 1
+
+    var J = JInput
+
+    J = if (settings.iteration % (2 * Npar) == 0 || dX2 > 0) {
+        lmFdJ(func, t, p, yHat, dp, settings).as2D() // finite difference
+    }
+    else {
+        lmBroydenJ(pOld, yOld, J, p, yHat).as2D() // rank-1 update
+    }
+
+    val deltaY = yDat.minus(yHat)
+
+    val chiSq = deltaY.transposed().dot( deltaY.times(weight) ).as2D()
+    val JtWJ = J.transposed().dot ( J.times( weight.dot(BroadcastDoubleTensorAlgebra.ones(ShapeND(intArrayOf(1, Npar)))) ) ).as2D()
+    val JtWdy = J.transposed().dot( weight.times(deltaY) ).as2D()
+
+    return arrayOf(JtWJ,JtWdy,chiSq,yHat,J)
+}
+
+private fun lmBroydenJ(pOld: MutableStructure2D<Double>, yOld: MutableStructure2D<Double>, JInput: MutableStructure2D<Double>,
+                       p: MutableStructure2D<Double>, y: MutableStructure2D<Double>): MutableStructure2D<Double> {
+    var J = JInput.copyToTensor()
+
+    val h = p.minus(pOld)
+    val increase = y.minus(yOld).minus( J.dot(h) ).dot(h.transposed()).div( (h.transposed().dot(h)).as2D()[0, 0] )
+    J = J.plus(increase)
+
+    return J.as2D()
+}
+
+private fun lmFdJ(func: (MutableStructure2D<Double>, MutableStructure2D<Double>, exampleNumber: Int) -> MutableStructure2D<Double>,
+                  t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, y: MutableStructure2D<Double>,
+                  dp: MutableStructure2D<Double>, settings: LMSettings): MutableStructure2D<Double> {
+    // default: dp = 0.001 * ones(1,n)
+
+    val m = length(y) // number of data points
+    val n = length(p) // number of parameters
+
+    val ps = p.copyToTensor().as2D()
+    val J = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, n))).as2D()  // initialize Jacobian to Zero
+    val del = BroadcastDoubleTensorAlgebra.zeros(ShapeND(intArrayOf(n, 1))).as2D()
+
+    for (j in 0 until n) {
+
+        del[j, 0] = dp[j, 0] * (1 + kotlin.math.abs(p[j, 0])) // parameter perturbation
+        p[j, 0] = ps[j, 0] + del[j, 0] // perturb parameter p(j)
+
+        val epsilon = 0.0000001
+        if (kotlin.math.abs(del[j, 0]) > epsilon) {
+            val y1 = evaluateFunction(func, t, p, settings.exampleNumber)
+            settings.funcCalls += 1
+
+            if (dp[j, 0] < 0) { // backwards difference
+                for (i in 0 until J.shape.component1()) {
+                    J[i, j] = (y1.as2D().minus(y).as2D())[i, 0] / del[j, 0]
+                }
+            }
+            else {
+                // Do tests for it
+                p[j, 0] = ps[j, 0] - del[j, 0] // central difference, additional func call
+                for (i in 0 until J.shape.component1()) {
+                    J[i, j] = (y1.as2D().minus(evaluateFunction(func, t, p, settings.exampleNumber)).as2D())[i, 0] / (2 * del[j, 0])
+                }
+                settings.funcCalls += 1
+            }
+        }
+
+        p[j, 0] = ps[j, 0]
+    }
+
+    return J.as2D()
+}
diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
index cf3697e76..086c69e49 100644
--- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
@@ -7,12 +7,10 @@ package space.kscience.kmath.tensors.core.internal
 
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.IntBuffer
-import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.indices
+import space.kscience.kmath.structures.*
 import space.kscience.kmath.tensors.core.*
 import kotlin.math.abs
+import kotlin.math.max
 import kotlin.math.min
 import kotlin.math.sqrt
 
@@ -308,3 +306,299 @@ internal fun DoubleTensorAlgebra.svdHelper(
         matrixV.source[i] = vBuffer[i]
     }
 }
+
+private fun pythag(a: Double, b: Double): Double {
+    val at: Double = abs(a)
+    val bt: Double = abs(b)
+    val ct: Double
+    val result: Double
+    if (at > bt) {
+        ct = bt / at
+        result = at * sqrt(1.0 + ct * ct)
+    } else if (bt > 0.0) {
+        ct = at / bt
+        result = bt * sqrt(1.0 + ct * ct)
+    } else result = 0.0
+    return result
+}
+
+private fun SIGN(a: Double, b: Double): Double {
+    if (b >= 0.0)
+        return abs(a)
+    else
+        return -abs(a)
+}
+
+internal fun MutableStructure2D<Double>.svdGolubKahanHelper(u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
+                                                            v: MutableStructure2D<Double>, iterations: Int, epsilon: Double) {
+    val shape = this.shape
+    val m = shape.component1()
+    val n = shape.component2()
+    var f = 0.0
+    val rv1 = DoubleArray(n)
+    var s = 0.0
+    var scale = 0.0
+    var anorm = 0.0
+    var g = 0.0
+    var l = 0
+
+    val wStart = 0
+    val wBuffer = w.source
+
+    for (i in 0 until n) {
+        /* left-hand reduction */
+        l = i + 1
+        rv1[i] = scale * g
+        g = 0.0
+        s = 0.0
+        scale = 0.0
+        if (i < m) {
+            for (k in i until m) {
+                scale += abs(this[k, i]);
+            }
+            if (abs(scale) > epsilon) {
+                for (k in i until m) {
+                    this[k, i] = (this[k, i] / scale)
+                    s += this[k, i] * this[k, i]
+                }
+                f = this[i, i]
+                if (f >= 0) {
+                    g = (-1) * abs(sqrt(s))
+                } else {
+                    g = abs(sqrt(s))
+                }
+                val h = f * g - s
+                this[i, i] = f - g
+                if (i != n - 1) {
+                    for (j in l until n) {
+                        s = 0.0
+                        for (k in i until m) {
+                            s += this[k, i] * this[k, j]
+                        }
+                        f = s / h
+                        for (k in i until m) {
+                            this[k, j] += f * this[k, i]
+                        }
+                    }
+                }
+                for (k in i until m) {
+                    this[k, i] = this[k, i] * scale
+                }
+            }
+        }
+
+        wBuffer[wStart + i] = scale * g
+        /* right-hand reduction */
+        g = 0.0
+        s = 0.0
+        scale = 0.0
+        if (i < m && i != n - 1) {
+            for (k in l until n) {
+                scale += abs(this[i, k])
+            }
+            if (abs(scale) > epsilon) {
+                for (k in l until n) {
+                    this[i, k] = this[i, k] / scale
+                    s += this[i, k] * this[i, k]
+                }
+                f = this[i, l]
+                if (f >= 0) {
+                    g = (-1) * abs(sqrt(s))
+                } else {
+                    g = abs(sqrt(s))
+                }
+                val h = f * g - s
+                this[i, l] = f - g
+                for (k in l until n) {
+                    rv1[k] = this[i, k] / h
+                }
+                if (i != m - 1) {
+                    for (j in l until m) {
+                        s = 0.0
+                        for (k in l until n) {
+                            s += this[j, k] * this[i, k]
+                        }
+                        for (k in l until n) {
+                            this[j, k] += s * rv1[k]
+                        }
+                    }
+                }
+                for (k in l until n) {
+                    this[i, k] = this[i, k] * scale
+                }
+            }
+        }
+        anorm = max(anorm, (abs(wBuffer[wStart + i]) + abs(rv1[i])));
+    }
+
+    for (i in n - 1 downTo 0) {
+        if (i < n - 1) {
+            if (abs(g) > epsilon) {
+                for (j in l until n) {
+                    v[j, i] = (this[i, j] / this[i, l]) / g
+                }
+                for (j in l until n) {
+                    s = 0.0
+                    for (k in l until n)
+                        s += this[i, k] * v[k, j]
+                    for (k in l until n)
+                        v[k, j] += s * v[k, i]
+                }
+            }
+            for (j in l until n) {
+                v[i, j] = 0.0
+                v[j, i] = 0.0
+            }
+        }
+        v[i, i] = 1.0
+        g = rv1[i]
+        l = i
+    }
+
+    for (i in min(n, m) - 1 downTo 0) {
+        l = i + 1
+        g =  wBuffer[wStart + i]
+        for (j in l until n) {
+            this[i, j] = 0.0
+        }
+        if (abs(g) > epsilon) {
+            g = 1.0 / g
+            for (j in l until n) {
+                s = 0.0
+                for (k in l until m) {
+                    s += this[k, i] * this[k, j]
+                }
+                f = (s / this[i, i]) * g
+                for (k in i until m) {
+                    this[k, j] += f * this[k, i]
+                }
+            }
+            for (j in i until m) {
+                this[j, i] *= g
+            }
+        } else {
+            for (j in i until m) {
+                this[j, i] = 0.0
+            }
+        }
+        this[i, i] += 1.0
+    }
+
+    var flag = 0
+    var nm = 0
+    var c = 0.0
+    var h = 0.0
+    var y = 0.0
+    var z = 0.0
+    var x = 0.0
+    for (k in n - 1 downTo 0) {
+        for (its in 1 until iterations) {
+            flag = 1
+            for (newl in k downTo 0) {
+                nm = newl - 1
+                if (abs(rv1[newl]) + anorm == anorm) {
+                    flag = 0
+                    l = newl
+                    break
+                }
+                if (abs(wBuffer[wStart + nm]) + anorm == anorm) {
+                    l = newl
+                    break
+                }
+            }
+
+            if (flag != 0) {
+                c = 0.0
+                s = 1.0
+                for (i in l until k + 1) {
+                    f = s * rv1[i]
+                    rv1[i] = c * rv1[i]
+                    if (abs(f) + anorm == anorm) {
+                        break
+                    }
+                    g = wBuffer[wStart + i]
+                    h = pythag(f, g)
+                    wBuffer[wStart + i] = h
+                    h = 1.0 / h
+                    c = g * h
+                    s = (-f) * h
+                    for (j in 0 until m) {
+                        y = this[j, nm]
+                        z = this[j, i]
+                        this[j, nm] = y * c + z * s
+                        this[j, i] = z * c - y * s
+                    }
+                }
+            }
+
+            z =  wBuffer[wStart + k]
+            if (l == k) {
+                if (z < 0.0) {
+                    wBuffer[wStart + k] = -z
+                    for (j in 0 until n)
+                        v[j, k] = -v[j, k]
+                }
+                break
+            }
+
+            x =  wBuffer[wStart + l]
+            nm = k - 1
+            y =  wBuffer[wStart + nm]
+            g = rv1[nm]
+            h = rv1[k]
+            f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
+            g = pythag(f, 1.0)
+            f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x
+            c = 1.0
+            s = 1.0
+
+            var i = 0
+            for (j in l until nm + 1) {
+                i = j + 1
+                g = rv1[i]
+                y = wBuffer[wStart + i]
+                h = s * g
+                g = c * g
+                z = pythag(f, h)
+                rv1[j] = z
+                c = f / z
+                s = h / z
+                f = x * c + g * s
+                g = g * c - x * s
+                h = y * s
+                y *= c
+
+                for (jj in 0 until n) {
+                    x = v[jj, j];
+                    z = v[jj, i];
+                    v[jj, j] = x * c + z * s;
+                    v[jj, i] = z * c - x * s;
+                }
+                z = pythag(f, h)
+                wBuffer[wStart + j] = z
+                if (abs(z) > epsilon) {
+                    z = 1.0 / z
+                    c = f * z
+                    s = h * z
+                }
+                f = c * g + s * y
+                x = c * y - s * g
+                for (jj in 0 until m) {
+                    y = this[jj, j]
+                    z = this[jj, i]
+                    this[jj, j] = y * c + z * s
+                    this[jj, i] = z * c - y * s
+                }
+            }
+            rv1[l] = 0.0
+            rv1[k] = f
+            wBuffer[wStart + k] = x
+        }
+    }
+
+    for (i in 0 until n) {
+        for (j in 0 until m) {
+            u[j, i] = this[j, i]
+        }
+    }
+}
\ No newline at end of file
diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt
index e5dc55f68..88ffb0bfe 100644
--- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt
+++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorOps.kt
@@ -212,6 +212,36 @@ public fun DoubleTensorAlgebra.svd(
     return Triple(uTensor.transposed(), sTensor, vTensor.transposed())
 }
 
+public fun DoubleTensorAlgebra.svdGolubKahan(
+    structureND: StructureND<Double>,
+    iterations: Int = 30, epsilon: Double = 1e-10
+): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
+    val size = structureND.dimension
+    val commonShape = structureND.shape.slice(0 until size - 2)
+    val (n, m) = structureND.shape.slice(size - 2 until size)
+    val uTensor = zeros(commonShape + intArrayOf(n, m))
+    val sTensor = zeros(commonShape + intArrayOf(m))
+    val vTensor = zeros(commonShape + intArrayOf(m, m))
+
+    val matrices = structureND.asDoubleTensor().matrices
+    val uTensors = uTensor.matrices
+    val sTensorVectors = sTensor.vectors
+    val vTensors = vTensor.matrices
+
+    for (index in matrices.indices) {
+        val matrix = matrices[index]
+        val matrixSize = matrix.shape.linearSize
+        val curMatrix = DoubleTensor(
+            matrix.shape,
+            matrix.source.view(0, matrixSize).copy()
+        )
+        curMatrix.as2D().svdGolubKahanHelper(uTensors[index].as2D(), sTensorVectors[index], vTensors[index].as2D(),
+            iterations, epsilon)
+    }
+
+    return Triple(uTensor, sTensor, vTensor)
+}
+
 /**
  * Returns eigenvalues and eigenvectors of a real symmetric matrix input or a batch of real symmetric matrices,
  * represented by a pair `eigenvalues to eigenvectors`.
diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt
index cae01bed8..7222fc7a6 100644
--- a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt
+++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleTensorAlgebra.kt
@@ -6,9 +6,7 @@
 package space.kscience.kmath.tensors.core
 
 
-import space.kscience.kmath.nd.ShapeND
-import space.kscience.kmath.nd.contentEquals
-import space.kscience.kmath.nd.get
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.testutils.assertBufferEquals
 import kotlin.test.Test
diff --git a/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt
new file mode 100644
index 000000000..4b031eb11
--- /dev/null
+++ b/kmath-tensors/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestLmAlgorithm.kt
@@ -0,0 +1,280 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.tensors.core
+
+import space.kscience.kmath.nd.MutableStructure2D
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.component1
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.max
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.plus
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.pow
+import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times
+import kotlin.math.roundToLong
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class TestLmAlgorithm {
+    companion object {
+        fun funcEasyForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+            val m = t.shape.component1()
+            var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))
+
+            if (exampleNumber == 1) {
+                yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0]))).times(p[0, 0]) + t.times(p[2, 0]).times(
+                    DoubleTensorAlgebra.exp((t.times(-1.0 / p[3, 0])))
+                )
+            }
+            else if (exampleNumber == 2) {
+                val mt = t.max()
+                yHat = (t.times(1.0 / mt)).times(p[0, 0]) +
+                        (t.times(1.0 / mt)).pow(2).times(p[1, 0]) +
+                        (t.times(1.0 / mt)).pow(3).times(p[2, 0]) +
+                        (t.times(1.0 / mt)).pow(4).times(p[3, 0])
+            }
+            else if (exampleNumber == 3) {
+                yHat = DoubleTensorAlgebra.exp((t.times(-1.0 / p[1, 0])))
+                    .times(p[0, 0]) + DoubleTensorAlgebra.sin((t.times(1.0 / p[3, 0]))).times(p[2, 0])
+            }
+
+            return yHat.as2D()
+        }
+
+        fun funcMiddleForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+            val m = t.shape.component1()
+            var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1)))
+
+            val mt = t.max()
+            for(i in 0 until p.shape.component1()){
+                yHat += (t.times(1.0 / mt)).times(p[i, 0])
+            }
+
+            for(i in 0 until 5){
+                yHat = funcEasyForLm(yHat.as2D(), p, exampleNumber).asDoubleTensor()
+            }
+
+            return yHat.as2D()
+        }
+
+        fun funcDifficultForLm(t: MutableStructure2D<Double>, p: MutableStructure2D<Double>, exampleNumber: Int): MutableStructure2D<Double> {
+            val m = t.shape.component1()
+            var yHat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf (m, 1)))
+
+            val mt = t.max()
+            for(i in 0 until p.shape.component1()){
+                yHat = yHat.plus( (t.times(1.0 / mt)).times(p[i, 0]) )
+            }
+
+            for(i in 0 until 4){
+                yHat = funcEasyForLm((yHat.as2D() + t).as2D(), p, exampleNumber).asDoubleTensor()
+            }
+
+            return yHat.as2D()
+        }
+    }
+    @Test
+    fun testLMEasy() = DoubleTensorAlgebra {
+        val lmMatxYDat = doubleArrayOf(
+            19.6594, 18.6096, 17.6792, 17.2747, 16.3065, 17.1458, 16.0467, 16.7023, 15.7809, 15.9807,
+            14.7620, 15.1128, 16.0973, 15.1934, 15.8636, 15.4763, 15.6860, 15.1895, 15.3495, 16.6054,
+            16.2247, 15.9854, 16.1421, 17.0960, 16.7769, 17.1997, 17.2767, 17.5882, 17.5378, 16.7894,
+            17.7648, 18.2512, 18.1581, 16.7037, 17.8475, 17.9081, 18.3067, 17.9632, 18.2817, 19.1427,
+            18.8130, 18.5658, 18.0056, 18.4607, 18.5918, 18.2544, 18.3731, 18.7511, 19.3181, 17.3066,
+            17.9632, 19.0513, 18.7528, 18.2928, 18.5967, 17.8567, 17.7859, 18.4016, 18.9423, 18.4959,
+            17.8000, 18.4251, 17.7829, 17.4645, 17.5221, 17.3517, 17.4637, 17.7563, 16.8471, 17.4558,
+            17.7447, 17.1487, 17.3183, 16.8312, 17.7551, 17.0942, 15.6093, 16.4163, 15.3755, 16.6725,
+            16.2332, 16.2316, 16.2236, 16.5361, 15.3721, 15.3347, 15.5815, 15.6319, 14.4538, 14.6044,
+            14.7665, 13.3718, 15.0587, 13.8320, 14.7873, 13.6824, 14.2579, 14.2154, 13.5818, 13.8157
+        )
+
+        var exampleNumber = 1
+        val p_init = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(4, 1)), doubleArrayOf(5.0, 2.0, 0.2, 10.0)
+        ).as2D()
+
+        var t = ones(ShapeND(intArrayOf(100, 1))).as2D()
+        for (i in 0 until 100) {
+            t[i, 0] = t[i, 0] * (i + 1)
+        }
+
+        val yDat = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(100, 1)), lmMatxYDat
+        ).as2D()
+
+        val weight = 4.0
+
+        val dp = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+        ).as2D()
+
+        val pMin = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(4, 1)), doubleArrayOf(-50.0, -20.0, -2.0, -100.0)
+        ).as2D()
+
+        val pMax = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(4, 1)), doubleArrayOf(50.0, 20.0, 2.0, 100.0)
+        ).as2D()
+
+        val inputData = LMInput(::funcEasyForLm, p_init, t, yDat, weight, dp, pMin, pMax, 100,
+            doubleArrayOf(1e-3, 1e-3, 1e-1, 1e-1), doubleArrayOf(1e-2, 11.0, 9.0), 1, 10, exampleNumber)
+
+        val result = levenbergMarquardt(inputData)
+        assertEquals(13, result.iterations)
+        assertEquals(31, result.funcCalls)
+        assertEquals(0.9131368192633, (result.resultChiSq * 1e13).roundToLong() / 1e13)
+        assertEquals(3.7790980 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13)
+        assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters)
+        val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(4, 1)), doubleArrayOf(20.527230909086, 9.833627103230, 0.997571256572, 50.174445822506)
+        ).as2D()
+        result.resultParameters = result.resultParameters.map { x -> (x * 1e12).toLong() / 1e12}.as2D()
+        val receivedParameters = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(4, 1)), doubleArrayOf(result.resultParameters[0, 0], result.resultParameters[1, 0],
+                result.resultParameters[2, 0], result.resultParameters[3, 0])
+        ).as2D()
+        assertEquals(expectedParameters[0, 0], receivedParameters[0, 0])
+        assertEquals(expectedParameters[1, 0], receivedParameters[1, 0])
+        assertEquals(expectedParameters[2, 0], receivedParameters[2, 0])
+        assertEquals(expectedParameters[3, 0], receivedParameters[3, 0])
+    }
+
+    @Test
+    fun TestLMMiddle() = DoubleTensorAlgebra {
+        val NData = 100
+        val tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+        for (i in 0 until NData) {
+            tExample[i, 0] = tExample[i, 0] * (i + 1)
+        }
+
+        val Nparams = 20
+        val pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            pExample[i, 0] = pExample[i, 0] + i - 25
+        }
+
+        val exampleNumber = 1
+
+        val yHat =  funcMiddleForLm(tExample, pExample, exampleNumber)
+
+        val pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            pInit[i, 0] = (pExample[i, 0] + 0.9)
+        }
+
+        val t = tExample
+        val yDat = yHat
+        val weight = 1.0
+        val dp = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+        ).as2D()
+        var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+        pMin = pMin.div(1.0 / -50.0)
+        val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+        pMin = pMin.div(1.0 / 50.0)
+        val opts = doubleArrayOf(3.0, 7000.0, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 11.0, 9.0, 1.0)
+
+        val inputData = LMInput(::funcMiddleForLm,
+            pInit.as2D(),
+            t,
+            yDat,
+            weight,
+            dp,
+            pMin.as2D(),
+            pMax.as2D(),
+            opts[1].toInt(),
+            doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+            doubleArrayOf(opts[6], opts[7], opts[8]),
+            opts[9].toInt(),
+            10,
+            1)
+
+        val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+        assertEquals(46, result.iterations)
+        assertEquals(113, result.funcCalls)
+        assertEquals(0.000005977, (result.resultChiSq * 1e9).roundToLong() / 1e9)
+        assertEquals(1.0 * 1e-7, (result.resultLambda * 1e13).roundToLong() / 1e13)
+        assertEquals(result.typeOfConvergence, TypeOfConvergence.InReducedChiSquare)
+        val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf( -23.9717, -18.6686, -21.7971,
+                -20.9681, -22.086, -20.5859, -19.0384, -17.4957, -15.9991, -14.576, -13.2441, -
+                12.0201, -10.9256, -9.9878, -9.2309, -8.6589, -8.2365, -7.8783, -7.4598, -6.8511)).as2D()
+        result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D()
+        val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            receivedParameters[i, 0] = result.resultParameters[i, 0]
+            assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0])
+        }
+    }
+
+    @Test
+    fun TestLMDifficult() = DoubleTensorAlgebra {
+        val NData = 200
+        var tExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(NData, 1))).as2D()
+        for (i in 0 until NData) {
+            tExample[i, 0] = tExample[i, 0] * (i + 1) - 104
+        }
+
+        val Nparams = 15
+        var pExample = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            pExample[i, 0] = pExample[i, 0] + i - 25
+        }
+
+        val exampleNumber = 1
+
+        var yHat =  funcDifficultForLm(tExample, pExample, exampleNumber)
+
+        var pInit = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            pInit[i, 0] = (pExample[i, 0] + 0.9)
+        }
+
+        var t = tExample
+        val yDat = yHat
+        val weight = 1.0 / Nparams * 1.0 - 0.085
+        val dp = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(1, 1)), DoubleArray(1) { -0.01 }
+        ).as2D()
+        var pMin = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+        pMin = pMin.div(1.0 / -50.0)
+        val pMax = DoubleTensorAlgebra.ones(ShapeND(intArrayOf(Nparams, 1)))
+        pMin = pMin.div(1.0 / 50.0)
+        val opts = doubleArrayOf(3.0, 7000.0, 1e-2, 1e-3, 1e-2, 1e-2, 1e-2, 11.0, 9.0, 1.0)
+
+        val inputData = LMInput(::funcDifficultForLm,
+            pInit.as2D(),
+            t,
+            yDat,
+            weight,
+            dp,
+            pMin.as2D(),
+            pMax.as2D(),
+            opts[1].toInt(),
+            doubleArrayOf(opts[2], opts[3], opts[4], opts[5]),
+            doubleArrayOf(opts[6], opts[7], opts[8]),
+            opts[9].toInt(),
+            10,
+            1)
+
+        val result = DoubleTensorAlgebra.levenbergMarquardt(inputData)
+
+        assertEquals(2375, result.iterations)
+        assertEquals(4858, result.funcCalls)
+        assertEquals(5.14347, (result.resultLambda * 1e5).roundToLong() / 1e5)
+        assertEquals(result.typeOfConvergence, TypeOfConvergence.InParameters)
+        val expectedParameters = BroadcastDoubleTensorAlgebra.fromArray(
+            ShapeND(intArrayOf(Nparams, 1)), doubleArrayOf(-23.6412, -16.7402, -21.5705, -21.0464,
+                -17.2852, -17.2959, -17.298, 0.9999, -17.2885, -17.3008, -17.2941, -17.2923, -17.2976, -17.3028, -17.2891)).as2D()
+        result.resultParameters = result.resultParameters.map { x -> (x * 1e4).roundToLong() / 1e4}.as2D()
+        val receivedParameters = zeros(ShapeND(intArrayOf(Nparams, 1))).as2D()
+        for (i in 0 until Nparams) {
+            receivedParameters[i, 0] = result.resultParameters[i, 0]
+            assertEquals(expectedParameters[i, 0], result.resultParameters[i, 0])
+        }
+    }
+}
\ No newline at end of file