\section{LEFTDIVIDE Matrix Equation Solver/Divide Operator}

\subsection{Usage}

The divide operator \verb|\| is really a combination of three
operators, all of which have the same general syntax:
\begin{verbatim}
  Y = A \ B
\end{verbatim}
where \verb|A| and \verb|B| are arrays of numerical type.  The result \verb|Y| depends
on which of the following three situations applies to the arguments
\verb|A| and \verb|B|:
\begin{enumerate}
\item  \verb|A| is a scalar, \verb|B| is an arbitrary \verb|n|-dimensional numerical array, in which case the output is each element of \verb|B| divided by the scalar \verb|A|.

\item  \verb|B| is a scalar, \verb|A| is an arbitrary \verb|n|-dimensional numerical array, in which case the output is the scalar \verb|B| divided by each element of \verb|A|.

\item  \verb|A,B| are matrices with the same number of rows, i.e., \verb|A| is of size \verb|M x K|, and \verb|B| is of size \verb|M x L|, in which case the output is of size \verb|K x L|.

\end{enumerate}
The output follows the standard type promotion rules, although in the first two cases, if \verb|A| and \verb|B| are integers, the output is an integer also, while in the third case if \verb|A| and \verb|B| are integers, the output is of type \verb|double|.

A few additional words about the third version, in which \verb|A| and \verb|B| are matrices.  Very loosely speaking, \verb|Y| is the matrix that satisfies \verb|A * Y = B|.  In cases where such a matrix exists.  If such a matrix does not exist, then a matrix \verb|Y| is returned that approximates \verb|A * Y \approx B|.
\subsection{Function Internals}

There are three formulae for the times operator.  For the first form
\[
Y(m_1,\ldots,m_d) = \frac{B(m_1,\ldots,m_d)}{A},
\]
and the second form
\[
Y(m_1,\ldots,m_d) = \frac{B}{A(m_1,\ldots,m_d)}.
\]
In the third form, the calculation of the output depends on the size of \verb|A|. Because each column of \verb|B| is treated independantly, we can rewrite the equation \verb|A Y = B| as
\[
  A [y_1, y_2,\ldots, y_l] = [b_1, b_2, \ldots, b_l]
\]
where \verb|y_i| are the columns of \verb|Y|, and \verb|b_i| are the columns of the matrix \verb|B|. If \verb|A| is a square matrix, then the LAPACK routine \verb|*gesvx| (where the \verb|*| is replaced with \verb|sdcz| depending on the type of the arguments) is used, which uses an LU decomposition of \verb|A| to solve the sequence of equations sequentially.  If \verb|A| is singular, then a warning is emitted. 

On the other hand, if \verb|A| is rectangular, then the LAPACK routine \verb|*gelsy| is used.  Note that these routines are designed to work with matrices \verb|A| that are full rank - either full column rank or full row rank.  If \verb|A| fails to satisfy this assumption, a warning is emitted.  If \verb|A| has full column rank (and thus necessarily has more rows than columns), then theoretically, this operator finds the columns \verb|y_i| that satisfy:
\[
  y_i = \arg \min_y \| A y - b_i \|_2
\]
and each column is thus the Least Squares solution of \verb|A y = b_i|.  On the other hand, if \verb|A| has full row rank (and thus necessarily has more columns than rows), then theoretically, this operator finds the columns \verb|y_i| that satisfy
\[
  y_i = \arg \min_{A y = b_i} \| y \|_2
\]
and each column is thus the Minimum Norm vector \verb|y_i| that satisfies \verb|A y_i = b_i|.  
In the event that the matrix \verb|A| is neither full row rank nor full column rank, a solution is returned, that is the minimum norm least squares solution.  The solution is computed using an orthogonal factorization technique that is documented in the LAPACK User's Guide (see the References section for details).
\subsection{Examples}

Here are some simple examples of the divide operator.  We start with a simple example of a full rank, square matrix:
\begin{verbatim}
--> A = [1,1;0,1]

A = 

 1 1 
 0 1 


\end{verbatim}
Suppose we wish to solve
\[
  \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}
  \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}
 = 
  \begin{bmatrix} 3 \\ 2 \end{bmatrix}
\]
(which by inspection has the solution \verb|y_1 = 1|, \verb|y_2 = 2|).  Thus we compute:
\begin{verbatim}
--> B = [3;2]

B = 

 3 
 2 

--> Y = A\B

Y = 

 1 
 2 


\end{verbatim}

Suppose we wish to solve a trivial Least Squares (LS) problem.  We want to find a simple scaling of the vector \verb|[1;1]| that is closest to the point \verb|[2,1]|.  This is equivalent to solving
\[
\begin{bmatrix} 1 \\ 1 \end{bmatrix} y = \begin{bmatrix} 2 \\ 1 \end{bmatrix}
\]
in a least squares sense.  For fun, we can calculate the solution using calculus by hand.  The error we wish to minimize is
\[
  \varepsilon(y) = (y - 2)^2 + (y-1)^2.
\]
Taking a derivative with respect to \verb|y|, and setting to zero (which we must have for an extrema when \verb|y| is unconstrained)
\[
  2 (y-2) + 2 (y-1) = 0
\]
which we can simplify to \verb|4y = 6| or \verb|y = 3/2| (we must, technically, check to make sure this is a minimum, and not a maximum or an inflection point).  Here is the same calculation performed using FreeMat:
\begin{verbatim}
--> A = [1;1]

A = 

 1 
 1 

--> B = [2;1]

B = 

 2 
 1 

--> A\B

ans = 

    1.5000 


\end{verbatim}
which is the same solution.
