added solve method. Refactored blockLU to return inverses to new solv…

…e method. Refactored shiftIndices method. Updated docs.

docs/mllib-guide.md

@@ -65,6 +65,7 @@ We list major functionality from both below, with links to detailed guides.
 * [PMML model export](mllib-pmml-model-export.html)
 * [Matrix Factorization](mllib-matrix-factorization.html)
   * [LU Factorization](mllib-matrix-factorization.html#lu-factorization)
+  * [Solving a Linear System of Equations](mllib-matrix-factorization.html#solving-a-linear-system-of-equations)
 
 MLlib is under active development.
 The APIs marked `Experimental`/`DeveloperApi` may change in future releases,

docs/mllib-matrix-factorization.md

@@ -9,8 +9,9 @@ displayTitle: <a href="mllib-guide.html">MLlib</a> - Matrix Factorization
 
 #LU Factorization
 
-[`blockLU()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) returns a [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix)
-The LU decomposition is a well studied algorithm for decomposing a matrix into a lower diagonal matrix $L$ and an upper diagonal matrix $U$.  
+
+The LU decomposition is a well studied algorithm for decomposing a matrix into a lower diagonal matrix 
+$L$ and an upper diagonal matrix $U$, given as  
 
 `\begin{equation}
     PA = LU\\
@@ -32,10 +33,19 @@ The LU decomposition is a well studied algorithm for decomposing a matrix into a
         \label{eq:generalLUFactorization}
 \end{equation}`
 
-and $P$ is a row permutation matrix.  This algorithm is a highly stable method for inverting a matrix and solving linear systems of equations that appear in machine learning applications.  Larger linear equations of the type $Ax=b$ are usually solved with SGD, BGFS, or other gradient based methods.  Being able to solve these equations at scale to numerical precision rather than to convergence should open up possiblities for new algorithms within MLlib as well as other applications.  It scales as ~$n^3$, for a [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) with $n \times n$ blocks.
+where $P$ is a row permutation matrix. [`blockLU()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) 
+returns a Tuple containing the $P,L$ and $U$ [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) 
+objects.
 
-Once the decomposition is computed, the inverse is straightforward to compute using a forward subsitution method on $L$ and $U$.  The inversion method is not yet available in [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix).
+This algorithm is a highly stable method for inverting a matrix and solving linear systems of equations that appear in
+machine learning applications.  Larger linear equations of the type $AX=B$ are usually solved with SGD, BGFS, or other
+gradient based methods.  Solving these equations at scale to numerical precision should open up possiblities for new algorithms within MLlib, as well as other applications.  
+It scales as ~$N^3$, for a [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) 
+with $N \times N$ blocks.
 
+Once the decomposition is computed, many other quantities can be easily derived.  The most important one,   
+[`BlockMatrix.solve()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix), and is 
+described in the next section.
 
 **Example**
 
@@ -73,7 +83,7 @@ println( "error (sb ~6.7e-16): " + error.toString )
 
 **How it Works**
 
-The LU decomposition of $A$ can be written in 4 block form as:
+The LU decomposition of $A$ can be written as four block matrices, given as
 `\begin{align}
     PA & = LU\\
         \begin{pmatrix}
@@ -91,7 +101,7 @@ The LU decomposition of $A$ can be written in 4 block form as:
         & =         \begin{pmatrix}
             L_{11}U_{11}&L_{11}U_{12} \\ 
             L_{21}U_{11}&L_{21}U_{12}+L_{22}U_{22}
-            \end{pmatrix}
+            \end{pmatrix}.
         \label{eq:basicLUBlockDecomposition}
 \end{align}`
 
@@ -104,26 +114,36 @@ P_2 A_{21} & = L_{21}U_{11}               & \Rightarrow  & L_{21}
 P_2 A_{22} & = L_{21}U_{12}+L_{22}U_{22}  & \Rightarrow  & (P_{2},L_{22},U_{22})  & = & \text{LU}_{RECURSIVE}(S) \label{eq:A22Solve}\\
 \end{align}`
 
-where $A_{11}$ is chosen to be a single block, so that the Breeze library can be called to compute $\eqref{eq:A11Solve}$.  The Breeze library will return the $P_{1},L_{11}$ and $U_{11}$ matrices.  Equation $\eqref{eq:A22Solve}$ is a recursive call that generates successive calculations of the Schur Complement, given as 
+where $A_{11}$ is chosen to be a single block, so that the Breeze library can be called to compute 
+$\eqref{eq:A11Solve}$.  The Breeze library will return the $P_{1},L_{11}$ and $U_{11}$ matrices.  
+
+
+Equation $\eqref{eq:A22Solve}$ is a recursive call that generates successive calculations of the Schur Complement, 
+given as 
 
 `\begin{equation}
 S = A_{22} - L_{21} U_{12}.
 \label{eq:SchurComplementInitialForm}
 \end{equation}`
 
-In this form, there is a dependency on the calculation of $\eqref{eq:A11Solve}$.  An equivalent form of the Schur complement can be used by substituting equations $\eqref{eq:U12Solve}$ and $\eqref{eq:L21Solve}$, giving
+In this form, there is a dependency on the calculation of $\eqref{eq:A11Solve}$.  An equivalent form of the Schur 
+complement can be used by substituting equations $\eqref{eq:U12Solve}$ and $\eqref{eq:L21Solve}$, giving
 
 `\begin{equation}
 S = A_{22} - A_{21} A_{11}^{-1} A_{12}, 
 \label{eq:SchurComplementFinalForm}
 \end{equation}`
 
-which allows for the calculation of equation $\eqref{eq:A22Solve}$ with no dependency on equation $\eqref{eq:A11Solve}$, resulting in a slight increase in parallelism at the expense of recomputing the inverse of $A_{11}$ on a separate process.  
-
-The Schur Complement in $\eqref{eq:A22Solve}$ is computed using $\eqref{eq:SchurComplementFinalForm}$ and passed recursively to the next itereration.  In this way, $P_2$ is never explicitly calculated, but is built up as a set of $P_1$ matrices.  The computation completes when the Schur Complement is a single block.
+which allows for the calculation of equation $\eqref{eq:A22Solve}$ with no dependency on equation 
+$\eqref{eq:A11Solve}$, resulting in a slight increase in parallelism at the expense of recomputing the inverse of
+$A_{11}$ on a separate process.  
 
-
-  Equations $\eqref{eq:U12Solve}$ and $\eqref{eq:L21Solve}$ are computed with [`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) operations.  The Schur Complement in $\eqref{eq:A22Solve}$ is computed using $\eqref{eq:SchurComplementFinalForm}$ with a [`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) opperation and passed recursively to the next itereration.  In this way, $P_2$ is never explicitly calculated, but is built up as a set of $P_1$ matrices.  The computation completes when the Schur Complement is a single block.
+Equations $\eqref{eq:U12Solve}$ and $\eqref{eq:L21Solve}$ are computed with 
+[`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) operations.  
+The Schur Complement in $\eqref{eq:A22Solve}$ is computed using $\eqref{eq:SchurComplementFinalForm}$ with a 
+[`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) 
+operation and passed recursively to the next iteration.  In this way, $P_2$ is never explicitly calculated, 
+but is built up as a set of $P_1$ matrices.  The computation completes when the Schur Complement is a single block.
 
 Instead of building the solution incrementally, and using more frequent but smaller block matrix operations, we construct large [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) structures, and carry out the multiplication at the end of the calculation.  This should leverage the optimizations present in the [`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) routine more effectively.  The matrices formed to carry out the operations are described int the in the figure below.
 
@@ -138,4 +158,91 @@ Instead of building the solution incrementally, and using more frequent but smal
 
 The Matrix multiply operations shown in the figure are generalizations of equations $\eqref{eq:A11Solve}$, $\eqref{eq:U12Solve}$, and $\eqref{eq:L21Solve}$.
 
+#Solving a Linear System of Equations
+[`BlockMatrix.solve()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) calls the LU Factorization
+method to solve the linear system of equations $AX=B$ for $X$.
+The method of forward subsitution solution is used.  This method is is described 
+[`here`](https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution), for elements of a matrix, and 
+and the generalization to blocks using [`BlockMatrix.solve()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix)
+is described briefly here.  
+
+For a Linear system of equations $AX=B$, where $A$ is an $N \times N$ matrix, and $X$ and $B$ are  
+[`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) objects of with $N$ 
+row blocks and $W$ column blocks.  $B$ can have any number of columns, including 1.  The number of row and columns 
+per block must be the same, however.  [`BlockMatrix.solve(B)`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix)  
+will return $X$ as a [`BlockMatrix`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix) object
+with the same dimensions of $B$.
+
+The forward substitution equations are straightforward do obtain in two steps.  The equation $AX=PB$  
+can be expressed as $LUX=\widehat{B}$.  First, the solution to $LY=\widehat{B}$ for $Y$ is obtained, followed by 
+solving $UX=Y$ for $X$.
+
+Expanding $LY=\widehat{B}$ as [`BlockMatrix.multiply()`](api/scala/index.html#org.apache.spark.mllib.linalg.distributed.BlockMatrix)
+operations gives
+
+`\begin{align}
+\widehat{B}_m &=& \sum_{i=1}^{m}L_{m,i}Y_m \\
+ &=& L_{m,m}Y_m + \sum_{i=1}^{m-1}L_{m,i}Y_m. \\
+\end{align}`
+
+The solution for the $m$th row of $Y$ (which contains $W$ column elements) in terms of previous rows with the following
+recurrence relation
+
+`\begin{equation}
+Y_m = L_{m,m}^{-1}\left [  \widehat{B}_m-\sum_{i=1}^{m}L_{m,i}Y_m\right ].\\
+\end{equation}`
+
+Once $Y$ is obtained, a similar method is applied to solve $UX=Y$.  Here we build the solution from the last row block, 
+so the indexing is reversed.  Letting $\widetilde{m}=N-m+1$ and expanding the multiplication terms gives
+
+
+`\begin{align}
+Y_{\widetilde{m}} &=& \sum_{i=\widetilde{m}}^{N}U_{\widetilde{m},i}X_{\widetilde{m}} \\
+ &=& U_{\widetilde{m},\widetilde{m}}X_{\widetilde{m}} + \sum_{i=\widetilde{m}+1}^{N}U_{\widetilde{m},i}X_{\widetilde{m}}. \\
+\end{align}`
+
+$X$ is obtained, starting with the bottom row block as $\widetilde{m}=N$ and incrementing to $\widetilde{m}=1$, using  
+
+`\begin{equation}
+X_\widetilde{m} = U_{\widetilde{m},\widetilde{m}}^{-1}\left[ Y_{\widetilde{m}}-\sum_{i=\widetilde{m}+1}^{N}U_{\widetilde{m},i}X_\widetilde{m}\right ].\\
+\end{equation}`
 
+**Example**
+{% highlight scala %}
+import breeze.linalg.{DenseMatrix => BDM}
+import org.apache.spark.mllib.linalg._
+import org.apache.spark.mllib.linalg.distributed.BlockMatrix
+
+    // square matrix, but not fully populated in edge blocks.
+val blocksForLU: Seq[((Int, Int), Matrix)] = Seq(
+  ((0, 0), new DenseMatrix(2, 2, Array(1.0, 2.0, 3.0, 2.0))),
+  ((0, 1), new DenseMatrix(2, 2, Array(2.0, 1.0, 3.0, 5.0))),
+  ((1, 0), new DenseMatrix(2, 2, Array(3.0, 2.0, 1.0, 1.0))),
+  ((1, 1), new DenseMatrix(2, 2, Array(1.0, 2.0, 0.0, 1.0))),
+  ((0, 2), new DenseMatrix(2, 1, Array(1.0, 1.0))),
+  ((1, 2), new DenseMatrix(2, 1, Array(1.0, 3.0))),
+  ((2, 0), new DenseMatrix(1, 2, Array(1.0, 0.0))),
+  ((2, 1), new DenseMatrix(1, 2, Array(1.0, 2.0))),
+  ((2, 2), new DenseMatrix(1, 1, Array(4.0))))
+
+val A = new BlockMatrix(sc.parallelize(blocksForLU), 2, 2)
+
+val B = new BlockMatrix(sc.parallelize(Seq( //5x7 B vector (7 columns of 5 row B)
+      ((0, 0), new DenseMatrix(2, 2, Array( 0, 1, -1, 2))),
+      ((0, 1), new DenseMatrix(2, 2, Array(-2, -1, -3, -2))),
+      ((0, 2), new DenseMatrix(2, 2, Array(-4, -3, -5, -4))),
+      ((0, 3), new DenseMatrix(2, 2, Array(-5, -6, -7, -7))),
+      ((1, 0), new DenseMatrix(2, 2, Array( 2, 3, 1, 2))),
+      ((1, 1), new DenseMatrix(2, 2, Array( 3, 1, -1, 4))),
+      ((1, 2), new DenseMatrix(2, 2, Array(-2, -1, -3, -2))),
+      ((1, 3), new DenseMatrix(2, 2, Array(2, 1, 3, 2))),
+      ((2, 0), new DenseMatrix(1, 2, Array( 0, 5))),
+      ((2, 1), new DenseMatrix(1, 2, Array( 1, 3))),
+      ((2, 2), new DenseMatrix(1, 2, Array( 0, 0))),
+      ((2, 3), new DenseMatrix(1, 2, Array( 1, 1))))), 2, 2)
+
+val X = A.solve(B)
+val residual = A.multiply(X).subtract(B)
+val error = residual.toLocalMatrix.toArray.reduce(_ + Math.abs(_))
+println("error" + error.toString)//sb 2.25e-14
+{% endhighlight %}