diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt
index 3d1625a50..d980c510f 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/LinearOpsTensorAlgebra.kt
@@ -14,7 +14,7 @@ public interface LinearOpsTensorAlgebra<T, TensorType : TensorStructure<T>, Inde
     public fun TensorType.cholesky(): TensorType
 
     //https://pytorch.org/docs/stable/linalg.html#torch.linalg.qr
-    public fun TensorType.qr(): TensorType
+    public fun TensorType.qr(): Pair<TensorType, TensorType>
 
     //https://pytorch.org/docs/stable/generated/torch.lu.html
     public fun TensorType.lu(): Pair<TensorType, IndexTensorType>
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt
index 9a4d13d2c..d99c35569 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt
@@ -1,5 +1,7 @@
 package space.kscience.kmath.tensors.core
 
+import space.kscience.kmath.nd.MutableStructure1D
+import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.structures.*
 import space.kscience.kmath.tensors.TensorStructure
 
@@ -16,6 +18,9 @@ public open class BufferedTensor<T>(
     public val numel: Int
         get() = linearStructure.size
 
+    internal constructor(tensor: BufferedTensor<T>) :
+            this(tensor.shape, tensor.buffer, tensor.bufferStart)
+
     override fun get(index: IntArray): T = buffer[bufferStart + linearStructure.offset(index)]
 
     override fun set(index: IntArray, value: T) {
@@ -30,6 +35,39 @@ public open class BufferedTensor<T>(
 
     override fun hashCode(): Int = 0
 
+    internal fun vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
+        val n = shape.size
+        val vectorOffset = shape[n - 1]
+        val vectorShape = intArrayOf(shape.last())
+        for (offset in 0 until numel step vectorOffset) {
+            val vector = BufferedTensor(vectorShape, buffer, offset)
+            yield(vector)
+        }
+    }
+
+    internal fun matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
+        check(shape.size >= 2) { "todo" }
+        val n = shape.size
+        val matrixOffset = shape[n - 1] * shape[n - 2]
+        val matrixShape = intArrayOf(shape[n - 2], shape[n - 1])
+        for (offset in 0 until numel step matrixOffset) {
+            val matrix = BufferedTensor(matrixShape, buffer, offset)
+            yield(matrix)
+        }
+    }
+
+    internal inline fun forEachVector(vectorAction: (BufferedTensor<T>) -> Unit): Unit {
+        for (vector in vectorSequence()) {
+            vectorAction(vector)
+        }
+    }
+
+    internal inline fun forEachMatrix(matrixAction: (BufferedTensor<T>) -> Unit): Unit {
+        for (matrix in matrixSequence()) {
+            matrixAction(matrix)
+        }
+    }
+
 }
 
 public class IntTensor internal constructor(
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt
index b81cf0364..6be50e710 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleLinearOpsTensorAlgebra.kt
@@ -1,8 +1,11 @@
 package space.kscience.kmath.tensors.core
 
-import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
+import space.kscience.kmath.nd.MutableStructure1D
+import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.as1D
 import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.tensors.LinearOpsTensorAlgebra
+import kotlin.math.sqrt
 
 public class DoubleLinearOpsTensorAlgebra :
     LinearOpsTensorAlgebra<Double, DoubleTensor, IntTensor>,
@@ -12,6 +15,48 @@ public class DoubleLinearOpsTensorAlgebra :
 
     override fun DoubleTensor.det(): DoubleTensor = detLU()
 
+    private inline fun luHelper(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>, m: Int) {
+        for (row in 0 until m) pivots[row] = row
+
+        for (i in 0 until m) {
+            var maxVal = -1.0
+            var maxInd = i
+
+            for (k in i until m) {
+                val absA = kotlin.math.abs(lu[k, i])
+                if (absA > maxVal) {
+                    maxVal = absA
+                    maxInd = k
+                }
+            }
+
+            //todo check singularity
+
+            if (maxInd != i) {
+
+                val j = pivots[i]
+                pivots[i] = pivots[maxInd]
+                pivots[maxInd] = j
+
+                for (k in 0 until m) {
+                    val tmp = lu[i, k]
+                    lu[i, k] = lu[maxInd, k]
+                    lu[maxInd, k] = tmp
+                }
+
+                pivots[m] += 1
+
+            }
+
+            for (j in i + 1 until m) {
+                lu[j, i] /= lu[i, i]
+                for (k in i + 1 until m) {
+                    lu[j, k] -= lu[j, i] * lu[i, k]
+                }
+            }
+        }
+    }
+
     override fun DoubleTensor.lu(): Pair<DoubleTensor, IntTensor> {
 
         checkSquareMatrix(shape)
@@ -35,6 +80,37 @@ public class DoubleLinearOpsTensorAlgebra :
 
     }
 
+    private inline fun pivInit(
+        p: MutableStructure2D<Double>,
+        pivot: MutableStructure1D<Int>,
+        n: Int
+    ) {
+        for (i in 0 until n) {
+            p[i, pivot[i]] = 1.0
+        }
+    }
+
+    private inline fun luPivotHelper(
+        l: MutableStructure2D<Double>,
+        u: MutableStructure2D<Double>,
+        lu: MutableStructure2D<Double>,
+        n: Int
+    ) {
+        for (i in 0 until n) {
+            for (j in 0 until n) {
+                if (i == j) {
+                    l[i, j] = 1.0
+                }
+                if (j < i) {
+                    l[i, j] = lu[i, j]
+                }
+                if (j >= i) {
+                    u[i, j] = lu[i, j]
+                }
+            }
+        }
+    }
+
     override fun luPivot(
         luTensor: DoubleTensor,
         pivotsTensor: IntTensor
@@ -42,7 +118,8 @@ public class DoubleLinearOpsTensorAlgebra :
         //todo checks
         checkSquareMatrix(luTensor.shape)
         check(
-            luTensor.shape.dropLast(1).toIntArray() contentEquals pivotsTensor.shape
+            luTensor.shape.dropLast(2).toIntArray() contentEquals pivotsTensor.shape.dropLast(1).toIntArray() ||
+                    luTensor.shape.last() == pivotsTensor.shape.last() - 1
         ) { "Bed shapes ((" } //todo rewrite
 
         val n = luTensor.shape.last()
@@ -63,6 +140,27 @@ public class DoubleLinearOpsTensorAlgebra :
 
     }
 
+    private inline fun choleskyHelper(
+        a: MutableStructure2D<Double>,
+        l: MutableStructure2D<Double>,
+        n: Int
+    ) {
+        for (i in 0 until n) {
+            for (j in 0 until i) {
+                var h = a[i, j]
+                for (k in 0 until j) {
+                    h -= l[i, k] * l[j, k]
+                }
+                l[i, j] = h / l[j, j]
+            }
+            var h = a[i, i]
+            for (j in 0 until i) {
+                h -= l[i, j] * l[i, j]
+            }
+            l[i, i] = sqrt(h)
+        }
+    }
+
     override fun DoubleTensor.cholesky(): DoubleTensor {
         // todo checks
         checkSquareMatrix(shape)
@@ -76,8 +174,52 @@ public class DoubleLinearOpsTensorAlgebra :
         return lTensor
     }
 
-    override fun DoubleTensor.qr(): DoubleTensor {
-        TODO("ANDREI")
+    private fun MutableStructure1D<Double>.dot(other: MutableStructure1D<Double>): Double {
+        var res = 0.0
+        for (i in 0 until size) {
+            res += this[i] * other[i]
+        }
+        return res
+    }
+
+    private fun MutableStructure1D<Double>.l2Norm(): Double {
+        return sqrt((0 until size).sumOf { this[it] * this[it] })
+    }
+
+    fun qrHelper(
+        matrix: MutableStructure2D<Double>,
+        q: MutableStructure2D<Double>,
+        r: MutableStructure2D<Double>
+    ) {
+        //todo check square
+        val n = matrix.colNum
+        for (j in 0 until n) {
+            val v = matrix.columns[j]
+            if (j > 0) {
+                for (i in 0 until j) {
+                    r[i, j] = q.columns[i].dot(matrix.columns[j])
+                    for (k in 0 until n) {
+                        v[k] = v[k] - r[i, j] * q.columns[i][k]
+                    }
+                }
+            }
+            r[j, j] = v.l2Norm()
+            for (i in 0 until n) {
+                q[i, j] = v[i] / r[j, j]
+            }
+        }
+    }
+
+    override fun DoubleTensor.qr(): Pair<DoubleTensor, DoubleTensor> {
+        checkSquareMatrix(shape)
+        val qTensor = zeroesLike()
+        val rTensor = zeroesLike()
+        val seq = matrixSequence().zip((qTensor.matrixSequence().zip(rTensor.matrixSequence())))
+        for ((matrix, qr) in seq) {
+            val (q, r) = qr
+            qrHelper(matrix.as2D(), q.as2D(), r.as2D())
+        }
+        return qTensor to rTensor
     }
 
     override fun DoubleTensor.svd(): Triple<DoubleTensor, DoubleTensor, DoubleTensor> {
@@ -88,6 +230,14 @@ public class DoubleLinearOpsTensorAlgebra :
         TODO("ANDREI")
     }
 
+    private fun luMatrixDet(luTensor: MutableStructure2D<Double>, pivotsTensor: MutableStructure1D<Int>): Double {
+        val lu = luTensor.as2D()
+        val pivots = pivotsTensor.as1D()
+        val m = lu.shape[0]
+        val sign = if ((pivots[m] - m) % 2 == 0) 1.0 else -1.0
+        return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
+    }
+
     public fun DoubleTensor.detLU(): DoubleTensor {
         val (luTensor, pivotsTensor) = lu()
         val n = shape.size
@@ -108,6 +258,33 @@ public class DoubleLinearOpsTensorAlgebra :
         return detTensor
     }
 
+    private fun luMatrixInv(
+        lu: MutableStructure2D<Double>,
+        pivots: MutableStructure1D<Int>,
+        invMatrix: MutableStructure2D<Double>
+    ) {
+        val m = lu.shape[0]
+
+        for (j in 0 until m) {
+            for (i in 0 until m) {
+                if (pivots[i] == j) {
+                    invMatrix[i, j] = 1.0
+                }
+
+                for (k in 0 until i) {
+                    invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
+                }
+            }
+
+            for (i in m - 1 downTo 0) {
+                for (k in i + 1 until m) {
+                    invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
+                }
+                invMatrix[i, j] /= lu[i, i]
+            }
+        }
+    }
+
     public fun DoubleTensor.invLU(): DoubleTensor {
         val (luTensor, pivotsTensor) = lu()
         val invTensor = luTensor.zeroesLike()
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
index c878fc58c..b418f3647 100644
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
@@ -1,7 +1,8 @@
 package space.kscience.kmath.tensors.core
 
-import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
+import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.tensors.TensorPartialDivisionAlgebra
 import kotlin.math.abs
 
 public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, DoubleTensor> {
@@ -229,6 +230,23 @@ public open class DoubleTensorAlgebra : TensorPartialDivisionAlgebra<Double, Dou
         return this.view(other.shape)
     }
 
+    private inline fun dotHelper(
+        a: MutableStructure2D<Double>,
+        b: MutableStructure2D<Double>,
+        res: MutableStructure2D<Double>,
+        l: Int, m: Int, n: Int
+    ) {
+        for (i in 0 until l) {
+            for (j in 0 until n) {
+                var curr = 0.0
+                for (k in 0 until m) {
+                    curr += a[i, k] * b[k, j]
+                }
+                res[i, j] = curr
+            }
+        }
+    }
+
     override fun DoubleTensor.dot(other: DoubleTensor): DoubleTensor {
         if (this.shape.size == 1 && other.shape.size == 1) {
             return DoubleTensor(intArrayOf(1), doubleArrayOf(this.times(other).buffer.array().sum()))
diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt
deleted file mode 100644
index e78cc33c7..000000000
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/tensors/core/linutils.kt
+++ /dev/null
@@ -1,188 +0,0 @@
-package space.kscience.kmath.tensors.core
-
-import space.kscience.kmath.nd.MutableStructure1D
-import space.kscience.kmath.nd.MutableStructure2D
-import space.kscience.kmath.nd.as1D
-import space.kscience.kmath.nd.as2D
-import kotlin.math.sqrt
-
-
-internal inline fun <T> BufferedTensor<T>.vectorSequence(): Sequence<BufferedTensor<T>> = sequence {
-    val n = shape.size
-    val vectorOffset = shape[n - 1]
-    val vectorShape = intArrayOf(shape.last())
-    for (offset in 0 until numel step vectorOffset) {
-        val vector = BufferedTensor(vectorShape, buffer, offset)
-        yield(vector)
-    }
-}
-
-internal inline fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = sequence {
-    check(shape.size >= 2) { "todo" }
-    val n = shape.size
-    val matrixOffset = shape[n - 1] * shape[n - 2]
-    val matrixShape = intArrayOf(shape[n - 2], shape[n - 1])
-    for (offset in 0 until numel step matrixOffset) {
-        val matrix = BufferedTensor(matrixShape, buffer, offset)
-        yield(matrix)
-    }
-}
-
-internal inline fun <T> BufferedTensor<T>.forEachVector(vectorAction: (BufferedTensor<T>) -> Unit): Unit {
-    for (vector in vectorSequence()) {
-        vectorAction(vector)
-    }
-}
-
-internal inline fun <T> BufferedTensor<T>.forEachMatrix(matrixAction: (BufferedTensor<T>) -> Unit): Unit {
-    for (matrix in matrixSequence()) {
-        matrixAction(matrix)
-    }
-}
-
-
-internal inline fun dotHelper(
-    a: MutableStructure2D<Double>,
-    b: MutableStructure2D<Double>,
-    res: MutableStructure2D<Double>,
-    l: Int, m: Int, n: Int
-) {
-    for (i in 0 until l) {
-        for (j in 0 until n) {
-            var curr = 0.0
-            for (k in 0 until m) {
-                curr += a[i, k] * b[k, j]
-            }
-            res[i, j] = curr
-        }
-    }
-}
-
-internal inline fun luHelper(lu: MutableStructure2D<Double>, pivots: MutableStructure1D<Int>, m: Int) {
-    for (row in 0 until m) pivots[row] = row
-
-    for (i in 0 until m) {
-        var maxVal = -1.0
-        var maxInd = i
-
-        for (k in i until m) {
-            val absA = kotlin.math.abs(lu[k, i])
-            if (absA > maxVal) {
-                maxVal = absA
-                maxInd = k
-            }
-        }
-
-        //todo check singularity
-
-        if (maxInd != i) {
-
-            val j = pivots[i]
-            pivots[i] = pivots[maxInd]
-            pivots[maxInd] = j
-
-            for (k in 0 until m) {
-                val tmp = lu[i, k]
-                lu[i, k] = lu[maxInd, k]
-                lu[maxInd, k] = tmp
-            }
-
-            pivots[m] += 1
-
-        }
-
-        for (j in i + 1 until m) {
-            lu[j, i] /= lu[i, i]
-            for (k in i + 1 until m) {
-                lu[j, k] -= lu[j, i] * lu[i, k]
-            }
-        }
-    }
-}
-
-internal inline fun pivInit(
-    p: MutableStructure2D<Double>,
-    pivot: MutableStructure1D<Int>,
-    n: Int
-) {
-    for (i in 0 until n) {
-        p[i, pivot[i]] = 1.0
-    }
-}
-
-internal inline fun luPivotHelper(
-    l: MutableStructure2D<Double>,
-    u: MutableStructure2D<Double>,
-    lu: MutableStructure2D<Double>,
-    n: Int
-) {
-    for (i in 0 until n) {
-        for (j in 0 until n) {
-            if (i == j) {
-                l[i, j] = 1.0
-            }
-            if (j < i) {
-                l[i, j] = lu[i, j]
-            }
-            if (j >= i) {
-                u[i, j] = lu[i, j]
-            }
-        }
-    }
-}
-
-internal inline fun choleskyHelper(
-    a: MutableStructure2D<Double>,
-    l: MutableStructure2D<Double>,
-    n: Int
-) {
-    for (i in 0 until n) {
-        for (j in 0 until i) {
-            var h = a[i, j]
-            for (k in 0 until j) {
-                h -= l[i, k] * l[j, k]
-            }
-            l[i, j] = h / l[j, j]
-        }
-        var h = a[i, i]
-        for (j in 0 until i) {
-            h -= l[i, j] * l[i, j]
-        }
-        l[i, i] = sqrt(h)
-    }
-}
-
-internal inline fun luMatrixDet(luTensor: MutableStructure2D<Double>, pivotsTensor: MutableStructure1D<Int>): Double {
-    val lu = luTensor.as2D()
-    val pivots = pivotsTensor.as1D()
-    val m = lu.shape[0]
-    val sign = if ((pivots[m] - m) % 2 == 0) 1.0 else -1.0
-    return (0 until m).asSequence().map { lu[it, it] }.fold(sign) { left, right -> left * right }
-}
-
-internal inline fun luMatrixInv(
-    lu: MutableStructure2D<Double>,
-    pivots: MutableStructure1D<Int>,
-    invMatrix: MutableStructure2D<Double>
-) {
-    val m = lu.shape[0]
-
-    for (j in 0 until m) {
-        for (i in 0 until m) {
-            if (pivots[i] == j) {
-                invMatrix[i, j] = 1.0
-            }
-
-            for (k in 0 until i) {
-                invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
-            }
-        }
-
-        for (i in m - 1 downTo 0) {
-            for (k in i + 1 until m) {
-                invMatrix[i, j] -= lu[i, k] * invMatrix[k, j]
-            }
-            invMatrix[i, j] /= lu[i, i]
-        }
-    }
-}
diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt
index d2d4ffb67..b79c54dd1 100644
--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/tensors/core/TestDoubleLinearOpsAlgebra.kt
@@ -76,4 +76,50 @@ class TestDoubleLinearOpsTensorAlgebra {
         val b = fromArray(intArrayOf(3), doubleArrayOf(5.5, 2.6, 6.4))
         assertEquals(a.dot(b).value(), 59.92)
     }
+
+    @Test
+    fun testQR() = DoubleLinearOpsTensorAlgebra {
+        val shape = intArrayOf(2, 2, 2)
+        val buffer = doubleArrayOf(
+            1.0, 3.0,
+            1.0, 2.0,
+            1.5, 1.0,
+            10.0, 2.0
+        )
+
+        val tensor = fromArray(shape, buffer)
+
+        val (q, r) = tensor.qr()
+
+        assertTrue { q.shape contentEquals shape }
+        assertTrue { r.shape contentEquals shape }
+
+        assertTrue { q.dot(r).buffer.array().epsEqual(buffer) }
+
+        //todo check orthogonality/upper triang.
+    }
+
+    @Test
+    fun testLU() = DoubleLinearOpsTensorAlgebra {
+        val shape = intArrayOf(2, 2, 2)
+        val buffer = doubleArrayOf(
+            1.0, 3.0,
+            1.0, 2.0,
+            1.5, 1.0,
+            10.0, 2.0
+        )
+        val tensor = fromArray(shape, buffer)
+
+        val (lu, pivots) = tensor.lu()
+
+        // todo check lu
+
+        val (p, l, u) = luPivot(lu, pivots)
+
+        assertTrue { p.shape contentEquals shape }
+        assertTrue { l.shape contentEquals shape }
+        assertTrue { u.shape contentEquals shape }
+
+        assertTrue { p.dot(tensor).buffer.array().epsEqual(l.dot(u).buffer.array()) }
+    }
 }