Today, we will be looking at solving a simple linear least squares regression problem with Arrayfire. Least squares regression is a commonly used form of linear regression. The often large sample sizes and underlying matrix operations used for fitting Linear Least Squares makes it an excellent candidate for a parallelized gpu. The specific problem we will be examining today is the classic polynomial fit to a scattered collection of points.
<h3><a id="Mathematical_Background_4"></a>Mathematical Background</h3>
Linear regression methods model relationships between a dependent variable (<strong>y</strong>) and one or more independent variables (<strong>X</strong>). Linear regression explores the conditional probability ( <strong>y</strong> | <strong>X</strong> ) Data are modeled using <em>linear</em> predictor functions. Given this linearity assumption between the dependent and independent variables, and (n) data points of the form:
<strong>{y<sub>i</sub>,&nbsp;&nbsp;&nbsp;&nbsp; x<sub>i1</sub>, … , x<sub>ip</sub>}</strong> ,
we can construct a model of the form:

<strong>y<sub>i</sub> = b<sub>1</sub>x<sub>i1</sub> + … + b<sub>p</sub>x<sub>ip</sub> + <em>e</em><sub>i</sub></strong>,&nbsp;&nbsp;&nbsp;&nbsp; where <strong>i = 1, …, n</strong>.

or, if we stack each point’s linear model atop each other, in vector form:

<strong>y = Xb + <em>e</em></strong>    where,

<strong>y</strong> - vector of dependent variables
<strong>X</strong> - “Design Matrix” of size <em>(n x p)</em> . This matrix contains the independent variables.
<strong>b</strong> - parameter vector the regression will be fitting
<strong>e</strong> - error vector/noise

The goal of least squares linear regression is to find a parameter vector <strong>b</strong> such that the square of the error is minimized:

<strong>b = argmin<sub>b</sub> S( b ) = argmin<sub>b</sub>|| e || <sup>2</sup> = argmin<sub>b</sub> || y - Xb || <sup>2</sup></strong>,

where <strong>S( b )</strong> – cost function.

Solving the minimization can be done through various gradient descent methods which iteratively work to set the gradient vector of the cost function to <strong>0</strong>. The minimized gradient of the cost function can also be found analytically. The resulting equations are called the normal equations and can be written in vector form as:

<strong>(X<sup>T</sup>X) b = X<sup>T</sup>y</strong>

After a quick re-arrangement, we can algebraically solve for our parameter vector:

<strong>b = (X<sup>T</sup>X)<sup>-1</sup> X<sup>T</sup>y</strong>

On first glance, our specific problem of fitting an n-order polynomial(for example n = 5) might not seem like a linear problem since we have quadratic, cubic, and higher terms: <strong>y = x<sup>5</sup> + x<sup>4</sup> + x<sup>3</sup> + x<sup>2</sup> + x<sup>1</sup> + x<sup>0</sup></strong>, however we are actually fitting the parameter vector <em>b</em> and our equation is of the linear combination form:
<strong>y = b<sub>1</sub>x<sup>5</sup> + b<sub>2</sub>x<sup>4</sup> + b<sub>3</sub>x<sup>3</sup> + b<sub>4</sub>x<sup>2</sup> + b<sub>5</sub>x<sup>1</sup> + b<sub>6</sub>x<sup>0</sup></strong>
<h3><a id="Coding_with_Arrayfire_32"></a>Coding with Arrayfire</h3>
The set of points to be fitted will be located within a window of size (WIDTH, HEIGHT). First, generate a set of points within the window:
<pre><code>array particles = randu(NSAMPLES, 2); //generate (NSAMPLE x 2) random uniform matrix
//scale [0,1] to [0, WIDTH] and [0, HEIGHT]
particles.col(0) *= WIDTH;
particles.col(1) *= HEIGHT;</code></pre>
Then, form our regressors matrix <strong>( X )</strong>:
<pre><code>array X = constant( 1, NSAMPLES, 2, f32);
X.col(1) = particles.col(0); // X [NSAMPLES 2] = [ 1  xvals ]</code></pre>
With the dependent and independent variables prepared, evaluate the algebraic solution for the parameter vector <strong>b</strong>:
<pre><code>array b = matmul(  matmul(inverse(matmul(X.T(), X)), X.T()), // (X'X) * X'                                
                   particles.col(1) );                    //   * y</code></pre>
Now that we have found <strong>b</strong>, we have a complete linear model that predicts <strong>( y | X )</strong>. It would be nice if we could see this linear model, so we will attempt to draw a simple grayscale visualization. First, set up an “image” matrix and visualize our original dataset:
<pre><code>array image = constant(0, HEIGHT, WIDTH, f32); //declare a black image to match our window
array ids(NSAMPLES, u32); //indexes specifying where in the image we should draw white
ids = (particles.col(0).as(u32) * (HEIGHT)) + particles.col(1).as(u32);
image(ids) = 1; //fill the specified indices with white

array mask_point = array(5, 5, point_kernel);
image = dilate(image, mask_point); //dilate operation with point_kernel to prettify points</code></pre>
Then plot the fitted line on our “image”:
<pre><code>float *params = b.host(); //copy params from device to host array
array yplot = range(dim4(WIDTH)) * params[1] + params[0]; //y = b1*x + b0
delete [] params;  //implicitly allocated params must be freed

ids = (range(dim4(WIDTH)).as(u32) * (HEIGHT)) + yplot.as(u32); //computes the indices in our image that the line passes through
ids = ids(yplot &gt; 0 &amp;&amp; yplot &lt; HEIGHT); //filters the ids where our y-values are out of bounds</code></pre>
And voila!
<h3> <a href="http://sandbox.arrayfire.com/wp-content/uploads/2015/11/Screenshot-from-2015-11-12-16_56_511.png"><img class=" size-medium wp-image-8644 aligncenter" src="http://sandbox.arrayfire.com/wp-content/uploads/2015/11/Screenshot-from-2015-11-12-16_56_511-300x150.png" alt="Least square linear regression" width="300" height="150" /></a></h3>
<h3><a id="Coding_for_higher_order_polynomials_72"></a>Coding for higher order polynomials</h3>
Ideally, we would be able to fit any n-th order polynomial to our data, so first, prepare our design matrix to follow the higher order terms:
<pre><code>X = constant(1, NSAMPLES, N_order, f32);
 X.col(1) *= particles.col(0);
 for(int i=2; i&lt;N_order; ++i)
    X.col(i) = X.col(i-1) * particles.col(0); //replace each column with one order higher column
// X [NSAMPLES N_order] = [ 1 x x^2 x^3 ... x^N_order]</code></pre>
Our equations remain the same <strong>( (X<sup>T</sup>X) b = X<sup>T</sup>y )</strong>, and while the identical algebraic solution for <strong>b</strong> would work it is quite inefficient to find the inverse of a (now) large normal equations Gramian matrix <strong>( X<sup>T</sup>X )</strong>. Instead, since the Gramian matrix has the nice property of being a positive semi-definite matrix, we can substitute it with the result of a Cholesky decomposition. The equation slightly changes <strong>( U<sup>T</sup>U b = X<sup>T</sup>y )</strong> and the procedure is as follows:
<ul>
	<li>Find the upper triangular matrix <strong>( U )</strong> from the Cholesky decomposition of <strong>( X<sup>T</sup>X)</strong></li>
	<li>Solve for vector z in <strong>U<sup>T</sup>z = X<sup>T</sup>y</strong></li>
	<li>Finally, solve for the parameter vector <strong>b</strong> in <strong>Ub = z</strong></li>
</ul>
Suprisingly, Arrayfire has convenient functions for each of these steps.
<pre><code>array U;
cholesky(U, matmul(X.T(), X), true); //cholesky on X'X, true parameter for upper triangular portion
array z = solve(U.T(), matmul(X.T(), particles.col(1)), AF_MAT_LOWER); //solve and specify the coefficients are lower triangular for optimized solution method
b = solve(U, z, AF_MAT_UPPER); //solve and specify the coefficients are upper triangular</code></pre>
And our result:

<img class=" size-medium wp-image-8645 aligncenter" src="http://sandbox.arrayfire.com/wp-content/uploads/2015/11/Screenshot-from-2015-11-12-16_59_251-300x150.png" alt="Screenshot from 2015-11-12 16_59_25" width="300" height="150" />

The complete source code for this project can be found <a href="https://github.com/syurkevi/leastsquares_arrayfire_example">here.</a>
<h3><a id="Conclusion_98"></a>Conclusion</h3>
We have implemented a least squares linear regression with Arrayfire and also visualized its results. Arrayfire provides many convenient functions which allowed us to choose the most robust method of solving our problem without sacrificing convenience and ease of coding.