kmath_tensors_examples.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "@file:Repository(\"https://repo.kotlin.link\")\n",
    "@file:DependsOn(\"space.kscience:kmath-tensors-jvm:0.3.0-dev-8\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import space.kscience.kmath.tensors.core.*\n",
    "import space.kscience.kmath.operations.invoke"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "val randSeed = 100500L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take `x` values from range 0 to 9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 0.0      ,  1.0      ,  2.0      ,  3.0      ,  4.0      ,  5.0      ,  6.0      ,  7.0      ,  8.0      ,  9.0      ]\n",
       ")"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val x = DoubleTensorAlgebra {\n",
    "    fromArray(\n",
    "        intArrayOf(10),\n",
    "        (0 until 10).toList().map { it.toDouble() }.toDoubleArray()\n",
    "    )\n",
    "} \n",
    "\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define `y` as dependent on `x` with noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 2.4858   ,  6.0337   ,  7.7462   ,  9.3974   ,  9.9228   ,  1.1071e+1,  1.3791e+1,  1.5351e+1,  1.6051e+1,  2.0488e+1]\n",
       ")"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val y = DoubleTensorAlgebra { 2.0 * x + (3.0 + x.randomNormalLike(randSeed) * 1.5) }\n",
    "\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stack them into single `dataset` tensor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 0.0      ,  2.4858   ],\n",
       "    [ 1.0      ,  6.0337   ],\n",
       "    [ 2.0      ,  7.7462   ],\n",
       "    [ 3.0      ,  9.3974   ],\n",
       "    [ 4.0      ,  9.9228   ],\n",
       "    [ 5.0      ,  1.1071e+1],\n",
       "    [ 6.0      ,  1.3791e+1],\n",
       "    [ 7.0      ,  1.5351e+1],\n",
       "    [ 8.0      ,  1.6051e+1],\n",
       "    [ 9.0      ,  2.0488e+1]]\n",
       ")"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val dataset = DoubleTensorAlgebra { stack(listOf(x, y)).transpose() }\n",
    "dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalize `x` and `y` and save means and stds for further recovery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "val xMean = DoubleTensorAlgebra { x.mean() }\n",
    "val yMean = DoubleTensorAlgebra { y.mean() }\n",
    "\n",
    "val xStd = DoubleTensorAlgebra { x.std() }\n",
    "val yStd = DoubleTensorAlgebra { y.std() }\n",
    "\n",
    "val xScaled = DoubleTensorAlgebra { (x - xMean) / xStd }\n",
    "val yScaled = DoubleTensorAlgebra { (y - yMean) / yStd }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "val mean = DoubleTensorAlgebra { \n",
    "    fromArray(\n",
    "        intArrayOf(2),\n",
    "        doubleArrayOf(xMean, yMean)\n",
    "    ) \n",
    "}\n",
    "\n",
    "val std = DoubleTensorAlgebra {\n",
    "    fromArray(\n",
    "        intArrayOf(2),\n",
    "        doubleArrayOf(xStd, yStd)\n",
    "    ) \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 4.5      ,  1.1234e+1]\n",
       ")"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 3.0277   ,  5.304    ]\n",
       ")"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute covaration matrix and its eigenvector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 1.0      ,  9.8423e-1],\n",
       "    [ 9.8423e-1,  1.0      ]]\n",
       ")"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val covMatrix = DoubleTensorAlgebra { cov(listOf(xScaled, yScaled)) }\n",
    "\n",
    "covMatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [-7.0711e-1, -7.0711e-1]\n",
       ")"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val (_, evecs) = DoubleTensorAlgebra { covMatrix.symEig() }\n",
    "val eigVec = DoubleTensorAlgebra { evecs[0] }\n",
    "\n",
    "eigVec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reduce dimension of `dataset`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 2.2172   ,  1.5107   ,  1.0488   ,  5.9515e-1,  2.9156e-1, -9.505e-2 , -6.912e-1 , -1.1327   , -1.4597   , -2.2848   ]\n",
       ")"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val datasetReduced = DoubleTensorAlgebra { eigVec dot stack(listOf(xScaled, yScaled)) }\n",
    "\n",
    "datasetReduced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 2)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.shape.joinToString(\", \", \"(\", \")\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datasetReduced.shape.joinToString(\", \", \"(\", \")\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Restore the 7th row of `dataset`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "val n = 7\n",
    "val restored = BroadcastDoubleTensorAlgebra { (datasetReduced[n] dot eigVec.view(intArrayOf(1, 2))) * std + mean }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 6.9251   ,  1.5482e+1]\n",
       ")"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "restored"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 7.0      ,  1.5351e+1]\n",
       ")"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DoubleTensorAlgebra { dataset[n] }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OLS with SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take `alpha` from normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 3.8081e-1,  3.7449   ,  4.7119   ,  5.6986   ,  9.4929   ]\n",
       ")"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val alpha = DoubleTensorAlgebra { \n",
    "    randomNormal(\n",
    "        intArrayOf(5),\n",
    "        randSeed\n",
    "    ) + fromArray(\n",
    "        intArrayOf(5),\n",
    "        doubleArrayOf(1.0, 2.5, 3.4, 5.0, 10.1)\n",
    "    ) \n",
    "}\n",
    "\n",
    "alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take sample of size 20 from normal distribution for `x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "val x = DoubleTensorAlgebra { \n",
    "    randomNormal(\n",
    "        intArrayOf(20, 5),\n",
    "        randSeed\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate `y` and add gaussian noise $N(0, 0.05)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "val y = DoubleTensorAlgebra { x dot alpha }\n",
    "DoubleTensorAlgebra { y += y.randomNormalLike(randSeed) * 0.05 }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute SVD of `x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "val (u, singValues, v) = DoubleTensorAlgebra { x.svd() }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have to make sure the singular values of the matrix are not close to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 6.9973   ,  4.8911   ,  4.2973   ,  3.5892   ,  2.5359   ]\n",
       ")"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "singValues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "${Σ}^{-1}$ matrix can be restored from singular values with `diagonalEmbedding()` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 1.4291e-1,  0.0      ,  0.0      ,  0.0      ,  0.0      ],\n",
       "    [ 0.0      ,  2.0445e-1,  0.0      ,  0.0      ,  0.0      ],\n",
       "    [ 0.0      ,  0.0      ,  2.327e-1 ,  0.0      ,  0.0      ],\n",
       "    [ 0.0      ,  0.0      ,  0.0      ,  2.7861e-1,  0.0      ],\n",
       "    [ 0.0      ,  0.0      ,  0.0      ,  0.0      ,  3.9433e-1]]\n",
       ")"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val sigmaInv = DoubleTensorAlgebra { diagonalEmbedding(singValues.map{ if (abs(it) < 1e-3) 0.0 else 1.0/it }) }\n",
    "\n",
    "sigmaInv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate `alpha` with OLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 3.9396e-1,  3.7266   ,  4.6911   ,  5.6899   ,  9.4909   ]\n",
       ")"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val alphaOLS = DoubleTensorAlgebra { v dot sigmaInv dot u.transpose() dot y }\n",
    "\n",
    "alphaOLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 3.8081e-1,  3.7449   ,  4.7119   ,  5.6986   ,  9.4929   ]\n",
       ")"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure out MSE of approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0310797370107518E-4"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fun mse(yTrue: DoubleTensor, yPred: DoubleTensor): Double = DoubleTensorAlgebra {\n",
    "    require(yTrue.shape.size == 1)\n",
    "    require(yTrue.shape contentEquals yPred.shape)\n",
    "\n",
    "    val diff = yTrue - yPred\n",
    "    ((diff dot diff) / yTrue.shape[0].toDouble()).value()\n",
    "}\n",
    "\n",
    "mse(alpha, alphaOLS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a linear system using the LU decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [-2.0      ,  1.5      ,  6.8      , -2.4      ]\n",
       ")"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val trueX = DoubleTensorAlgebra { \n",
    "    fromArray(\n",
    "        intArrayOf(4),\n",
    "        doubleArrayOf(-2.0, 1.5, 6.8, -2.4)\n",
    "    ) \n",
    "}\n",
    "\n",
    "trueX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 5.0e-1   ,  1.05e+1  ,  4.5      ,  1.0      ],\n",
       "    [ 8.5      ,  9.0e-1   ,  1.28e+1  ,  1.0e-1   ],\n",
       "    [ 5.56     ,  9.19     ,  7.62     ,  5.45     ],\n",
       "    [ 1.0      ,  2.0      , -3.0      , -2.5      ]]\n",
       ")"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val a = DoubleTensorAlgebra { \n",
    "    fromArray(\n",
    "        intArrayOf(4, 4),\n",
    "        doubleArrayOf(\n",
    "            0.5, 10.5, 4.5, 1.0,\n",
    "            8.5, 0.9, 12.8, 0.1,\n",
    "            5.56, 9.19, 7.62, 5.45,\n",
    "            1.0, 2.0, -3.0, -2.5\n",
    "        )\n",
    "    )\n",
    "}\n",
    "\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [ 4.295e+1 ,  7.115e+1 ,  4.1401e+1, -1.34e+1  ]\n",
       ")"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val b = DoubleTensorAlgebra { a dot trueX }\n",
    "\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate $P$, $L$ and $U$, such that $PA = LU$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "val (p, l, u) = DoubleTensorAlgebra { a.lu() }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 0.0      ,  1.0      ,  0.0      ,  0.0      ],\n",
       "    [ 1.0      ,  0.0      ,  0.0      ,  0.0      ],\n",
       "    [ 0.0      ,  0.0      ,  0.0      ,  1.0      ],\n",
       "    [ 0.0      ,  0.0      ,  1.0      ,  0.0      ]]\n",
       ")"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 1.0      ,  0.0      ,  0.0      ,  0.0      ],\n",
       "    [ 5.8824e-2,  1.0      ,  0.0      ,  0.0      ],\n",
       "    [ 1.1765e-1,  1.8131e-1,  1.0      ,  0.0      ],\n",
       "    [ 6.5412e-1,  8.2332e-1,  7.4013e-1,  1.0      ]]\n",
       ")"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 8.5      ,  9.0e-1   ,  1.28e+1  ,  1.0e-1   ],\n",
       "    [ 0.0      ,  1.0447e+1,  3.7471   ,  9.9412e-1],\n",
       "    [ 0.0      ,  0.0      , -5.1852   , -2.692    ],\n",
       "    [ 0.0      ,  0.0      ,  0.0      ,  6.5585   ]]\n",
       ")"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 8.5      ,  9.0e-1   ,  1.28e+1  ,  1.0e-1   ],\n",
       "    [ 5.0e-1   ,  1.05e+1  ,  4.5      ,  1.0      ],\n",
       "    [ 1.0      ,  2.0      , -3.0      , -2.5      ],\n",
       "    [ 5.56     ,  9.19     ,  7.62     ,  5.45     ]]\n",
       ")"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DoubleTensorAlgebra { p dot a }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [[ 8.5      ,  9.0e-1   ,  1.28e+1  ,  1.0e-1   ],\n",
       "    [ 5.0e-1   ,  1.05e+1  ,  4.5      ,  1.0      ],\n",
       "    [ 1.0      ,  2.0      , -3.0      , -2.5      ],\n",
       "    [ 5.56     ,  9.19     ,  7.62     ,  5.45     ]]\n",
       ")"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DoubleTensorAlgebra { l dot u }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DoubleTensorAlgebra { (p dot a) eq (l dot u) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Ax = b;$\n",
    "\n",
    "$PAx = Pb;$\n",
    "\n",
    "$LUx = Pb;$\n",
    "\n",
    "let $y = Ux$, then\n",
    "\n",
    "$Ly = Pb$ -- this system can be easily solved, since the matrix L is lower triangular;\n",
    "\n",
    "$Ux = y$ can be solved the same way, since the matrix U is upper triangular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function returns solution `x` of a system `lx = b`, `l` should be lower triangular."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "fun solveLT(l: DoubleTensor, b: DoubleTensor): DoubleTensor = DoubleTensorAlgebra {\n",
    "    val n = l.shape[0]\n",
    "    val x = zeros(intArrayOf(n))\n",
    "    for (i in 0 until n) {\n",
    "        x[intArrayOf(i)] = (b[intArrayOf(i)] -  l[i].dot(x).value()) / l[intArrayOf(i, i)]\n",
    "    }\n",
    "    x\n",
    "}\n",
    "\n",
    "val y = DoubleTensorAlgebra { solveLT(l, p dot b) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create `solveUT` for solving system with upper triangular matrix with using permutation matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "val revMat = DoubleTensorAlgebra { u.zeroesLike() }\n",
    "val n = revMat.shape[0]\n",
    "for (i in 0 until n) {\n",
    "    revMat[intArrayOf(i, n - 1 - i)] = 1.0\n",
    "}\n",
    "\n",
    "fun solveUT(u: DoubleTensor, b: DoubleTensor): DoubleTensor = DoubleTensorAlgebra {\n",
    "    revMat dot solveLT(\n",
    "        revMat dot u dot revMat, revMat dot b\n",
    "    ) \n",
    "}\n",
    "\n",
    "val x = solveUT(u, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [-2.0      ,  1.5      ,  6.8      , -2.4      ]\n",
       ")"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DoubleTensor(\n",
       "   [-2.0      ,  1.5      ,  6.8      , -2.4      ]\n",
       ")"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trueX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Kotlin",
   "language": "kotlin",
   "name": "kotlin"
  },
  "language_info": {
   "codemirror_mode": "text/x-kotlin",
   "file_extension": ".kt",
   "mimetype": "text/x-kotlin",
   "name": "kotlin",
   "nbconvert_exporter": "",
   "pygments_lexer": "kotlin",
   "version": "1.5.20-dev-3998"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}