Matrices and Fibonacci 1

When I was beginning to learn Scala, I decided one of my projects would be rewrite this stuff in Scala. The idea behind that is to compute Fibonacci Numbers using matrices. I won’t bother describing how it is done, as this is well covered on both of those links. Suffice to say that I need to computer the nth power of a 2x2 matrix of integers. This project was all but forgotten for a long time…

So it happened that just yesterday I was cleaning my to-do list, and I came upon that. I took a second look into it, and decided to take the plunge. Now, Raganwald uses an optimization which enables him to do the operations on vectors, but, by now, I actually thought the whole exercise was trivial and pointless.

I took a brief look into matrix libraries for Scala. I found the Scalala (as in Scala Linear Algebra) library, which does everything I need, very efficiently at that. My problem would be reduced to something like this (untested):

def fib(n) = n match {
  case 0 => throw new IllegalArgumentException
  case 1 => 1
  case _ =>
    val m = DenseMatrix(2,2)(1); m(1,1) = 0
    val mPowered = m :^ n
    mPowered(0,1)
}

Boring, wasn’t it? Yeah, I thought so too. Instead, I decided to write my own Matrix library. Pretty trivial stuff, but at least I’d take something out of it, I thought. So, I’ll talk about that: writing my Matrix library. If that will bore you, you might as well stop right here. I’ll show, though, some of the design decisions and constrains I faced. For someone new to Scala, it might be interesting read.

Well, then, where to start? With a class, of course. Ideally, I’d like to write something like this:

class Matrix[T]

This way, I could use matrices of Double for most use cases, of Ints or Longs for things like the Fibonacci problem, and BigDecimal for tough stuff. Unfortunately, we can’t do that. We need access to operators such as +, - and * over type T, and that means I need an upper bound for T defining these operators. Unfortunately, these types do not share an interface defining those elements. You see, they are not really “classes” as far as Java is concerned, and not being a proper class, they can’t share an abstract class or interface with these methods. Well, BigDecimal is a class, alright, but that’s it.

Because of this, I decided, instead, to define my class as MatrixInt. This was enough for Fibonacci, but later I created a MatrixDouble as well, just to finish implementing the most basic methods such a class would need.

There is a way around that. I can get away with a structural type, but that would be so slow as to defeat the whole purpose of a Matrix class. So I’ll skip it altogether.

The next design decision is how to store my data, and what constructor to use. Generally speaking, there are two strategies for storing a matrix: as a dense matrix or as a sparse matrix. A dense matrix is simply stored as an addressable, fixed-size, contiguously allocated buffer. Or, in other words, an array. A sparse matrix is used for matrices that are very large in their dimensions, but are mostly populated with zeros. Performance can go either way, depending on what kind of data and algorithms you are using those matrices for. For Fibonacci, dense is the way to go. Since it is also the easier way to go, I chose that way.

We also need to know how many dimensions the matrices may have. For simplicity sake, I chose bi-dimensional matrices, and that’s it. With all that in mind, my class thus began:

class MatrixInt(array: Array[Array[Int]]) {

I decided to take an array of proper dimensions as my constructor, instead of providing something special myself. With Scala 2.8, you can pass Array.ofDim(2,2) and be done with it. We’ll see something similar further down too.

Now, one important thing to consider is that an Array is a mutable object. I have decided to make my matrices immutable, so that I could use them just like I use Ints or Doubles. But, to ensure that is true, I can’t really on the array passed with the constructor, as it might change. So, first thing I do is copy it. Here:

  private def clone(a: Array[Array[Int]]) =
    for((i: Int) <- (0 until a.size).toArray)
      yield for((j: Int) <- (0 until a(0).size).toArray)
        yield a(i)(j)

  private val self = clone(array)

This clone array method will help in places other than the constructor. For instance, if the user wants the Array representation of the matrix, I have to clone my internal representation before passing it on. But there are more interesting things about it. Let’s consider it…

One thing to notice about it, is that instead of having a single for/yield statement, I have two, the second being fed to the first’s yield. It is built this way because I want to create an Array of Arrays. Each yield will result in one collection. Two collections, two yield statements.

Then we have ranges: (0 until a.size) and (0 until a(0).size). A statement like (0 until 5) will produce a range which, when iterated over, will yield 0, 1, 2, 3 and 4. It can do a few other things too, but let’s not be bothered by them. Note that if you want an “inclusive”, range, one which will include the upper limit (5, in the example), you use “to” instead of “until”. Or you create a new range from the existing one with “inclusive”.

Next, “.toArray”. Though a range will do as far as generating indices, I want yield to return an Array. If you want yield to return collection X, pass a collection X to the for statement. Therefore, “.toArray” to convert the range into an array.

The final thing to consider is the last yield statement. Yield a(i)(j). As I started writing my methods, I got very similar patterns, with only the last yield changing. So I started to consider how to reuse that code. The thing to consider, here, is that the intent of this double for/yield pattern is to create an array of arrays of the proper size, and populate it with a function of my choosing. Thinking about it that way, a refactoring might be done as follow:

  protected def createArrays(rows: Int, cols: Int, f: (Int, Int) => Int) =
    for((i: Int) <- (0 until rows).toArray)
      yield for((j: Int) <- (0 until cols).toArray)
        yield f(i,j)

  protected def clone(a: Array[Array[Int]]) =
    createArrays(a.size, a(0).size, (row, col) => a(row)(col))

Here I changed from private to protected, as I decided if I ever subclassed it, there would be no point in hiding these methods. I do not publish them, though, because they are helper methods, instead of methods acting upon the object itself.

Next, a.size and a(0).size were a bit cryptic, weren’t they? When we work on the matrices algorithms, it would be preferable to call them by what they are:

  val rows = array.size
  val cols = array(0).size

We would also like to be able to get a particular row or a particular column. It seems getting a particular row would be easy. Just return the subarray:

  def row(n : Int): Array[Int] = self(n)

There is a serious problem with that, though. I’m returning a mutable object from my own internal data! We could solve that with clone, but clone only worked for bi-dimensional arrays. Well, before we create new clone, let’s consider a simpler alternative:

  private def _row(n: Int): Array[Int] = self(n)
  def row(n: Int): Array[Int] = _row(n) map (x => x)

The private definition is here so we can avoid creating new arrays unnecessarily. Internally, we’ll make sure we never change the value of our arrays. The public method creates a new array by the magic of “map”. The method “map” is, in fact, used when you do a for/yield, along with flatMap when the for/yield has multiple iterators. So, in fact, the createArrays method is actually implemented as using map.

A map will take a collection and produce a new collection with a given mapping. In this case, our mapping is (x => x), or, in other words, given x, return x, producing a copy of that array. We define cols in a likewise fashion:

  private def _col(n: Int): Array[Int] = self map (_(n))
  def col(n: Int): Array[Int] = _col(n)

In this case, the private definition and the public definition are the same. Still, I decided to keep a private definition so that I could use _row and _col throughout my class. We’ll also define a method to access a particular element of our matrix:

  def apply(i: Int, j: Int): Int =
    if (i >= rows || j >= cols || i < 0 || j < 0)
      throw new IndexOutOfBoundsException
    else
      self(i)(j)

The method “apply” has a special meaning in Scala, as the compiler translates object(parms) into object.apply(parms). This means that I will be able to access a particular element of a matrix m with “m(i,j)”. Of course, because a matrix has bounds, I must test for them, and throw an appropriate exception if need be.

Note that I do not declare an @throws annotation, though I might if I wished to. Also, note that I’m returning the result of the if/else expression, and its type should be Int. But in the if statement, I’m throwing an exception, so what gives?

As it happens, throwing an exception has type Nothing, which is subtype to everything. I’m not sure if the compiler actually assigns a type to “throw”, as this is a keyword. But if I used Scala’s library “error” function instead, the question of the type would be important. For methods and functions which never return, such as error, declaring its type to be Nothing makes it possible to make use of them inside statements whose return type is important.

We’ll stop today with just three more definitions. Let’s multiply a matrix by an integer, producing a new matrix, as well as add and subtract two arrays. Here is how:

  def *(n: Int): MatrixInt = 
    new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) * n))

  def +(other: MatrixInt): MatrixInt =
    if (rows != other.rows || cols != other.cols)
      throw new IllegalArgumentException
    else
      new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) + other(row,col)))

  def -(other: MatrixInt): MatrixInt =
    if (rows != other.rows || cols != other.cols)
      throw new IllegalArgumentException
    else
      new MatrixInt(createArrays(rows, cols, (row, col) => this(row,col) - other(row, col)))