You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@systemml.apache.org by re...@apache.org on 2017/08/19 05:24:23 UTC
systemml git commit: [MINOR] DML Tips and Tricks Jupyter notebook
Repository: systemml
Updated Branches:
refs/heads/master 1adfc7266 -> 8e06ff3cc
[MINOR] DML Tips and Tricks Jupyter notebook
Project: http://git-wip-us.apache.org/repos/asf/systemml/repo
Commit: http://git-wip-us.apache.org/repos/asf/systemml/commit/8e06ff3c
Tree: http://git-wip-us.apache.org/repos/asf/systemml/tree/8e06ff3c
Diff: http://git-wip-us.apache.org/repos/asf/systemml/diff/8e06ff3c
Branch: refs/heads/master
Commit: 8e06ff3cca487715ae47da25cd896f59a7fcd817
Parents: 1adfc72
Author: Berthold Reinwald <re...@us.ibm.com>
Authored: Fri Aug 18 21:49:47 2017 -0700
Committer: Berthold Reinwald <re...@us.ibm.com>
Committed: Fri Aug 18 22:12:32 2017 -0700
----------------------------------------------------------------------
...DML Tips and Tricks (aka Fun With DML).ipynb | 753 +++++++++++++++++++
1 file changed, 753 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/systemml/blob/8e06ff3c/samples/jupyter-notebooks/DML Tips and Tricks (aka Fun With DML).ipynb
----------------------------------------------------------------------
diff --git a/samples/jupyter-notebooks/DML Tips and Tricks (aka Fun With DML).ipynb b/samples/jupyter-notebooks/DML Tips and Tricks (aka Fun With DML).ipynb
new file mode 100644
index 0000000..23d975a
--- /dev/null
+++ b/samples/jupyter-notebooks/DML Tips and Tricks (aka Fun With DML).ipynb
@@ -0,0 +1,753 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "1. [Value-based join of two Matrices](#JoinMatrices)\n",
+ "* [Filter Matrix to include only Frequent Column Values](#FilterMatrix)\n",
+ "* [Construct (sparse) Matrix from (rowIndex, colIndex, values) triplets](#Construct_sparse_Matrix)\n",
+ "* [Find and remove duplicates in columns or rows](#Find_and_remove_duplicates)\n",
+ "* [Set based Indexing](#Set_based_Indexing)\n",
+ "* [Group by Aggregate using Linear Algebra](#Multi_column_Sorting)\n",
+ "* [Cumulative Summation with Decay Multiplier](#CumSum_Product)\n",
+ "* [Invert Lower Triangular Matrix](#Invert_Lower_Triangular_Matrix)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "from systemml import MLContext, dml, jvm_stdout\n",
+ "ml = MLContext(sc)\n",
+ "\n",
+ "print (ml.buildTime())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Value-based join of two Matrices<a id=\"JoinMatrices\"/>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given matrix M1 and M2, join M1 on column 2 with M2 on column 2, and return matching rows of M1."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 55,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "M1 \n",
+ "1.000 1.000\n",
+ "2.000 3.000\n",
+ "3.000 3.000\n",
+ "4.000 4.000\n",
+ "5.000 3.000\n",
+ "6.000 4.000\n",
+ "7.000 1.000\n",
+ "8.000 2.000\n",
+ "9.000 1.000\n",
+ "\n",
+ "M2 \n",
+ "1.000 1.000\n",
+ "2.000 8.000\n",
+ "3.000 3.000\n",
+ "4.000 3.000\n",
+ "5.000 1.000\n",
+ "\n",
+ "M1[,2] joined with M2[,2], and return matching M1 rows\n",
+ "1.000 1.000\n",
+ "2.000 3.000\n",
+ "3.000 3.000\n",
+ "5.000 3.000\n",
+ "7.000 1.000\n",
+ "9.000 1.000\n",
+ "\n",
+ "SystemML Statistics:\n",
+ "Total execution time:\t\t0.001 sec.\n",
+ "Number of executed Spark inst:\t0.\n",
+ "\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "prog = \"\"\"\n",
+ "M1 = matrix ('1 1 2 3 3 3 4 4 5 3 6 4 7 1 8 2 9 1', rows = 9, cols = 2)\n",
+ "M2 = matrix ('1 1 2 8 3 3 4 3 5 1', rows = 5, cols = 2)\n",
+ "\n",
+ "I = rowSums (outer (M1[,2], t(M2[,2]), \"==\")) # I : indicator matrix for M1\n",
+ "M12 = removeEmpty (target = M1, margin = \"rows\", select = I) # apply filter to retrieve join result\n",
+ "\n",
+ "print (\"M1 \\n\" + toString(M1))\n",
+ "print (\"M2 \\n\" + toString(M2))\n",
+ "print (\"M1[,2] joined with M2[,2], and return matching M1 rows\\n\" + toString(M12))\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Filter Matrix to include only Frequent Column Values <a id=\"FilterMatrix\"/>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given a matrix, filter the matrix to only include rows with column values that appear more often than MinFreq."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 53,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "1.000 1.000\n",
+ "2.000 3.000\n",
+ "3.000 3.000\n",
+ "4.000 4.000\n",
+ "5.000 3.000\n",
+ "6.000 4.000\n",
+ "7.000 1.000\n",
+ "8.000 2.000\n",
+ "9.000 1.000\n",
+ "\n",
+ "1.000 1.000\n",
+ "2.000 3.000\n",
+ "3.000 3.000\n",
+ "5.000 3.000\n",
+ "7.000 1.000\n",
+ "9.000 1.000\n",
+ "\n",
+ "SystemML Statistics:\n",
+ "Total execution time:\t\t0.001 sec.\n",
+ "Number of executed Spark inst:\t0.\n",
+ "\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "prog = \"\"\"\n",
+ "MinFreq = 3 # minimum frequency of tokens\n",
+ "\n",
+ "M = matrix ('1 1 2 3 3 3 4 4 5 3 6 4 7 1 8 2 9 1', rows = 9, cols = 2)\n",
+ "\n",
+ "gM = aggregate (target = M[,2], groups = M[,2], fn = \"count\") # gM: group by and count (grouped matrix)\n",
+ "gv = cbind (seq(1,nrow(gM)), gM) # gv: add group values to counts (group values)\n",
+ "fg = removeEmpty (target = gv * (gv[,2] >= MinFreq), margin = \"rows\") # fg: filtered groups\n",
+ "I = rowSums (outer (M[,2] ,t(fg[,1]), \"==\")) # I : indicator of size M with filtered groups\n",
+ "fM = removeEmpty (target = M, margin = \"rows\", select = I) # FM: filter matrix\n",
+ "\n",
+ "print (toString(M))\n",
+ "print (toString(fM))\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Construct (sparse) Matrix from (rowIndex, colIndex, values) triplets<a id=\"Construct_sparse_Matrix\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given rowIndex, colIndex, and values as column vectors, construct (sparse) matrix."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "I = matrix (\"1 3 3 4 5\", rows = 5, cols = 1)\n",
+ "J = matrix (\"2 3 4 1 6\", rows = 5, cols = 1)\n",
+ "V = matrix (\"10 20 30 40 50\", rows = 5, cols = 1)\n",
+ "\n",
+ "M = table (I, J, V)\n",
+ "print (toString (M))\n",
+ "\"\"\"\n",
+ "ml.execute(dml(prog).output('M')).get('M').toNumPy()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Find and remove duplicates in columns or rows<a id=\"Find_and_remove_duplicates\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Assuming values are sorted."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "X = matrix (\"1 2 3 3 3 4 5 10\", rows = 8, cols = 1)\n",
+ "\n",
+ "I = rbind (matrix (1,1,1), (X[1:nrow (X)-1,] != X[2:nrow (X),])); # compare current with next value\n",
+ "res = removeEmpty (target = X, margin = \"rows\", select = I); # select where different\n",
+ "\"\"\"\n",
+ "ml.execute(dml(prog).output('res')).get('res').toNumPy()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### No assumptions on values."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "X = matrix (\"3 2 1 3 3 4 5 10\", rows = 8, cols = 1)\n",
+ "\n",
+ "I = aggregate (target = X, groups = X[,1], fn = \"count\") # group and count duplicates\n",
+ "res = removeEmpty (target = seq (1, max (X[,1])), margin = \"rows\", select = (I != 0)); # select groups\n",
+ "\"\"\"\n",
+ "ml.execute(dml(prog).output('res')).get('res').toNumPy()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Order the values and then remove duplicates."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "X = matrix (\"3 2 1 3 3 4 5 10\", rows = 8, cols = 1)\n",
+ "\n",
+ "X = order (target = X, by = 1) # order values\n",
+ "I = rbind (matrix (1,1,1), (X[1:nrow (X)-1,] != X[2:nrow (X),]));\n",
+ "res = removeEmpty (target = X, margin = \"rows\", select = I);\n",
+ "\"\"\"\n",
+ "ml.execute(dml(prog).output('res')).get('res').toNumPy()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Set based Indexing<a id=\"Set_based_Indexing\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given a matrix X, and a indicator matrix J with indices into X. \n",
+ "Use J to perform operation on X, e.g. add value 10 to cells in X indicated by J. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "X = matrix (1, rows = 1, cols = 100)\n",
+ "J = matrix (\"10 20 25 26 28 31 50 67 79\", rows = 1, cols = 9)\n",
+ "\n",
+ "res = X + table (matrix (1, rows = 1, cols = ncol (J)), J, 10)\n",
+ "\n",
+ "print (toString (res))\n",
+ "\"\"\"\n",
+ "ml.execute(dml(prog).output('res')).get('res').toNumPy()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "## Group by Aggregate using Linear Algebra<a id=\"Multi_column_Sorting\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given a matrix PCV as (Position, Category, Value), sort PCV by category, and within each category by value in descending order. Create indicator vector for category changes, create distinct categories, and perform linear algebra operations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = \"\"\"\n",
+ "C = matrix ('50 40 20 10 30 20 40 20 30', rows = 9, cols = 1) # category data\n",
+ "V = matrix ('20 11 49 33 94 29 48 74 57', rows = 9, cols = 1) # value data\n",
+ "\n",
+ "PCV = cbind (cbind (seq (1, nrow (C), 1), C), V); # PCV representation\n",
+ "PCV = order (target = PCV, by = 3, decreasing = TRUE, index.return = FALSE);\n",
+ "PCV = order (target = PCV, by = 2, decreasing = FALSE, index.return = FALSE);\n",
+ "\n",
+ "# Find all rows of PCV where the category has a new value, in comparison to the previous row\n",
+ "\n",
+ "is_new_C = matrix (1, rows = 1, cols = 1);\n",
+ "if (nrow (C) > 1) {\n",
+ " is_new_C = rbind (is_new_C, (PCV [1:nrow(C) - 1, 2] < PCV [2:nrow(C), 2]));\n",
+ "}\n",
+ "\n",
+ "# Associate each category with its index\n",
+ "\n",
+ "index_C = cumsum (is_new_C); # cumsum\n",
+ "\n",
+ "# For each category, compute:\n",
+ "# - the list of distinct categories\n",
+ "# - the maximum value for each category\n",
+ "# - 0-1 aggregation matrix that adds records of the same category\n",
+ "\n",
+ "distinct_C = removeEmpty (target = PCV [, 2], margin = \"rows\", select = is_new_C);\n",
+ "max_V_per_C = removeEmpty (target = PCV [, 3], margin = \"rows\", select = is_new_C);\n",
+ "C_indicator = table (index_C, PCV [, 1], max (index_C), nrow (C)); # table\n",
+ "\n",
+ "sum_V_per_C = C_indicator %*% V\n",
+ "\"\"\"\n",
+ "\n",
+ "res = ml.execute(dml(prog).output('PCV','distinct_C', 'max_V_per_C', 'C_indicator', 'sum_V_per_C'))\n",
+ "print (res.get('PCV').toNumPy())\n",
+ "print (res.get('distinct_C').toNumPy())\n",
+ "print (res.get('max_V_per_C').toNumPy())\n",
+ "print (res.get('C_indicator').toNumPy())\n",
+ "print (res.get('sum_V_per_C').toNumPy())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Cumulative Summation with Decay Multiplier<a id=\"CumSum_Product\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Given matrix X, compute:\n",
+ "\n",
+ " Y[i] = X[i]\n",
+ " + X[i-1] * C[i]\n",
+ " + X[i-2] * C[i] * C[i-1]\n",
+ " + X[i-3] * C[i] * C[i-1] * C[i-2]\n",
+ " + ...\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "cumsum_prod_def = \"\"\"\n",
+ "cumsum_prod = function (Matrix[double] X, Matrix[double] C, double start)\n",
+ " return (Matrix[double] Y)\n",
+ "# Computes the following recurrence in log-number of steps:\n",
+ "# Y [1, ] = X [1, ] + C [1, ] * start;\n",
+ "# Y [i+1, ] = X [i+1, ] + C [i+1, ] * Y [i, ]\n",
+ "{\n",
+ " Y = X; P = C; m = nrow(X); k = 1;\n",
+ " Y [1,] = Y [1,] + C [1,] * start;\n",
+ " while (k < m) {\n",
+ " Y [k + 1:m,] = Y [k + 1:m,] + Y [1:m - k,] * P [k + 1:m,];\n",
+ " P [k + 1:m,] = P [1:m - k,] * P [k + 1:m,];\n",
+ " k = 2 * k;\n",
+ " }\n",
+ "}\n",
+ "\"\"\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this example we use cumsum_prod for cumulative summation with \"breaks\", that is, multiple cumulative summations in one."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = cumsum_prod_def + \"\"\"\n",
+ "X = matrix (\"1 2 3 4 5 6 7 8 9\", rows = 9, cols = 1);\n",
+ "\n",
+ "#Zeros in C cause \"breaks\" that restart the cumulative summation from 0\n",
+ "C = matrix (\"0 1 1 0 1 1 1 0 1\", rows = 9, cols = 1);\n",
+ "\n",
+ "Y = cumsum_prod (X, C, 0);\n",
+ "\n",
+ "print (toString(Y))\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this example, we copy selected rows downward to all consecutive non-selected rows."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = cumsum_prod_def + \"\"\"\n",
+ "X = matrix (\"1 2 3 4 5 6 7 8 9\", rows = 9, cols = 1);\n",
+ "\n",
+ "# Ones in S represent selected rows to be copied, zeros represent non-selected rows\n",
+ "S = matrix (\"1 0 0 1 0 0 0 1 0\", rows = 9, cols = 1);\n",
+ "\n",
+ "Y = cumsum_prod (X * S, 1 - S, 0);\n",
+ "\n",
+ "print (toString(Y))\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This is a naive implementation of cumulative summation with decay multiplier."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "cumsum_prod_naive_def = \"\"\"\n",
+ "cumsum_prod_naive = function (Matrix[double] X, Matrix[double] C, double start)\n",
+ " return (Matrix[double] Y)\n",
+ "{\n",
+ " Y = matrix (0, rows = nrow(X), cols = ncol(X));\n",
+ " Y [1,] = X [1,] + C [1,] * start;\n",
+ " for (i in 2:nrow(X))\n",
+ " {\n",
+ " Y [i,] = X [i,] + C [i,] * Y [i - 1,]\n",
+ " }\n",
+ "}\n",
+ "\"\"\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "There is a significant performance difference between the <b>naive</b> implementation and the <b>tricky</b> implementation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = cumsum_prod_def + cumsum_prod_naive_def + \"\"\"\n",
+ "X = rand (rows = 20000, cols = 10, min = 0, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "C = rand (rows = 20000, cols = 10, min = 0, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "\n",
+ "Y1 = cumsum_prod_naive (X, C, 0.123);\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = cumsum_prod_def + cumsum_prod_naive_def + \"\"\"\n",
+ "X = rand (rows = 20000, cols = 10, min = 0, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "C = rand (rows = 20000, cols = 10, min = 0, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "\n",
+ "Y2 = cumsum_prod (X, C, 0.123);\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Invert Lower Triangular Matrix<a id=\"Invert_Lower_Triangular_Matrix\"></a>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this example, we invert a lower triangular matrix using a the following divide-and-conquer approach. Given lower triangular matrix L, we compute its inverse X which is also lower triangular by splitting both matrices in the middle into 4 blocks (in a 2x2 fashion), and multiplying them together to get the identity matrix:\n",
+ "\n",
+ "\\begin{equation}\n",
+ "L \\text{ %*% } X = \\left(\\begin{matrix} L_1 & 0 \\\\ L_2 & L_3 \\end{matrix}\\right)\n",
+ "\\text{ %*% } \\left(\\begin{matrix} X_1 & 0 \\\\ X_2 & X_3 \\end{matrix}\\right)\n",
+ "= \\left(\\begin{matrix} L_1 X_1 & 0 \\\\ L_2 X_1 + L_3 X_2 & L_3 X_3 \\end{matrix}\\right)\n",
+ "= \\left(\\begin{matrix} I & 0 \\\\ 0 & I \\end{matrix}\\right)\n",
+ "\\nonumber\n",
+ "\\end{equation}\n",
+ "\n",
+ "If we multiply blockwise, we get three equations: \n",
+ "\n",
+ "$\n",
+ "\\begin{equation}\n",
+ "L1 \\text{ %*% } X1 = 1\\\\ \n",
+ "L3 \\text{ %*% } X3 = 1\\\\\n",
+ "L2 \\text{ %*% } X1 + L3 \\text{ %*% } X2 = 0\\\\\n",
+ "\\end{equation}\n",
+ "$\n",
+ "\n",
+ "Solving these equation gives the following formulas for X:\n",
+ "\n",
+ "$\n",
+ "\\begin{equation}\n",
+ "X1 = inv(L1) \\\\\n",
+ "X3 = inv(L3) \\\\\n",
+ "X2 = - X3 \\text{ %*% } L2 \\text{ %*% } X1 \\\\\n",
+ "\\end{equation}\n",
+ "$\n",
+ "\n",
+ "If we already recursively inverted L1 and L3, we can invert L2. This suggests an algorithm that starts at the diagonal and iterates away from the diagonal, involving bigger and bigger blocks (of size 1, 2, 4, 8, etc.) There is a logarithmic number of steps, and inside each step, the inversions can be performed in parallel using a parfor-loop.\n",
+ "\n",
+ "Function \"invert_lower_triangular\" occurs within more general inverse operations and matrix decompositions. The divide-and-conquer idea allows to derive more efficient algorithms for other matrix decompositions.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "invert_lower_triangular_def = \"\"\"\n",
+ "invert_lower_triangular = function (Matrix[double] LI)\n",
+ " return (Matrix[double] LO)\n",
+ "{\n",
+ " n = nrow (LI);\n",
+ " LO = matrix (0, rows = n, cols = n);\n",
+ " LO = LO + diag (1 / diag (LI));\n",
+ " \n",
+ " k = 1;\n",
+ " while (k < n)\n",
+ " {\n",
+ " LPF = matrix (0, rows = n, cols = n);\n",
+ " parfor (p in 0:((n - 1) / (2 * k)), check = 0)\n",
+ " {\n",
+ " i = 2 * k * p;\n",
+ " j = i + k;\n",
+ " q = min (n, j + k);\n",
+ " if (j + 1 <= q) {\n",
+ " L1 = LO [i + 1:j, i + 1:j];\n",
+ " L2 = LI [j + 1:q, i + 1:j];\n",
+ " L3 = LO [j + 1:q, j + 1:q];\n",
+ " LPF [j + 1:q, i + 1:j] = -L3 %*% L2 %*% L1;\n",
+ " }\n",
+ " }\n",
+ " LO = LO + LPF;\n",
+ " k = 2 * k;\n",
+ " }\n",
+ "}\n",
+ "\"\"\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = invert_lower_triangular_def + \"\"\"\n",
+ "n = 1000;\n",
+ "A = rand (rows = n, cols = n, min = -1, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "Mask = cumsum (diag (matrix (1, rows = n, cols = 1)));\n",
+ "\n",
+ "L = (A %*% t(A)) * Mask; # Generate L for stability of the inverse\n",
+ "\n",
+ "X = invert_lower_triangular (L);\n",
+ "\n",
+ "print (\"Maximum difference between X %*% L and Identity = \" + max (abs (X %*% L - diag (matrix (1, rows = n, cols = 1)))));\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This is a naive implementation of inverting a lower triangular matrix."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "invert_lower_triangular_naive_def = \"\"\"\n",
+ "invert_lower_triangular_naive = function (Matrix[double] LI)\n",
+ " return (Matrix[double] LO)\n",
+ "{\n",
+ " n = nrow (LI);\n",
+ " LO = diag (matrix (1, rows = n, cols = 1));\n",
+ " for (i in 1:n - 1)\n",
+ " {\n",
+ " LO [i,] = LO [i,] / LI [i, i];\n",
+ " LO [i + 1:n,] = LO [i + 1:n,] - LI [i + 1:n, i] %*% LO [i,];\n",
+ " }\n",
+ " LO [n,] = LO [n,] / LI [n, n];\n",
+ "}\n",
+ "\"\"\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The naive implementation is significantly slower than the divide-and-conquer implementation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "prog = invert_lower_triangular_naive_def + \"\"\"\n",
+ "n = 1000;\n",
+ "A = rand (rows = n, cols = n, min = -1, max = 1, pdf = \"uniform\", sparsity = 1.0);\n",
+ "Mask = cumsum (diag (matrix (1, rows = n, cols = 1)));\n",
+ "\n",
+ "L = (A %*% t(A)) * Mask; # Generate L for stability of the inverse\n",
+ "\n",
+ "X = invert_lower_triangular_naive (L);\n",
+ "\n",
+ "print (\"Maximum difference between X %*% L and Identity = \" + max (abs (X %*% L - diag (matrix (1, rows = n, cols = 1)))));\n",
+ "\"\"\"\n",
+ "with jvm_stdout():\n",
+ " ml.execute(dml(prog))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 2",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.11"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}