Matrices and Fibonacci 3

Continuing our series on creating a Matrix class for Scala, now that we have the basics laid down, let’s go back to the Fibonacci problem. We need to do exponentiation of matrices. To do that, we first need to get multiplication going.

While I have no intention of explaining how to multiply matrices – the web is full of resources if you aren’t familiar with it – it’s important to restate the problem:

To find element of (row i, column j) of the matrix resulting from AxB, I need to multiply the nth element of row i of matrix A by the nth element of column j of matrix B, for all elements of said row and column, and then add the results. An initial, misguided, attempt might go like this:

(for(row <- 0 to rows
     col <- 0 to B.cols
    ) yield A(i,col)*B(row,j)) reduceLeft (_+_)

The problem here is that you are multiplying each element of the row by every element of the column. What we want to do is to “pair” the first element of the row with the first element of the column, the second element of the row with the second of the column, and so forth, and then multiply each pair. This pairing is called “zip” in Scala. You’ll find “zip” – and the occasional variation, such as zipWithIndex – on sequence-like collections, such as List and Array. The zipping thing goes like this:

_row(i) zip other._col(j)

That will create an Array of tuples of Int: Array[Tuple2[Int,Int]], or simply Array[(Int,Int)]. Tuples can be thought of as a small “package” of objects. It’s not a collection – you do not traverse it or iterate through it. Also, each component of it has its own type, as you can see. The only thing you can do with tuples, aside from passing it around, is accessing a tuple’s element. You do this with the methods “_1”, “_2”, “_3”, etc.

Next, we need to do _1*_2 for each element of that array. You can’t do "map (_._1 * _._2)", as this translates into "map ((x1, x2) => x1._1 * x2._2)", but map only passes a single argument: the tuple. You can do this, then: "map (t => t._1 * t._2)". You can, but I won’t… Instead, I’ll pass "(_*_)", an anonymous function that takes two parameters and return the result of their multiplication, to a factory of functions.

Factories, as anyone familiar with OOP might know, are methods which create new instances of objects. So, what that has to do with a function? Well, a function in Scala is an object that implements the method “apply”. You might wonder if MatrixInt isn’t a function, then, as it implements apply itself. In fact, it could be one, if we chose to extend MatrixInt with the trait Function2 (two arguments to apply).

At any rate, we’ll use a function factory to produce a new instance of a function, based on (_*_). This factory, called Function.tupled (method “tupled” on singleton object Function), takes a function which receives multiple parameters, and converts it into a function which receives a single tuple instead. After that, we can apply the reduceLeft. Our computation, then, will go like this:

this._row(i) zip other._col(j) map Function.tupled(_*_) reduceLeft (_+_)

The method, then, is thus defined:

  def *(other: MatrixInt): MatrixInt =
    if (cols != other.rows)
      throw new IllegalArgumentException
    else
      new MatrixInt(createArrays(rows, other.cols, (row, col) => {
        this._row(row) zip other._col(col) map Function.tupled(_*_) reduceLeft (_+_)
      }))

If we do away with the check to see if the matrices are compatible in size, and make allowance for the verbosity of the identifiers being used (like row and col instead of I and j), it comes down to a one-liner. Truly:

def *(o: MatrixInt) = new MatrixInt(createArrays(rows, o.cols, (i,j) => _row(i) zip o._col(j) map Function.tupled(_*_) reduceLeft (_+_)))

The only thing left to do is exponentiation. We could do it as a one-liner too:

def ^(n: Int) = (0 until n).foldLeft(unitMatrixInt(rows))((acc,_) => acc*this)

The method foldLeft is similar to reduce. But where reduce can be thought of as (((a+b)+c)+d)+e, foldLeft is more like this:

var acc = ?
acc = op(acc,a)
acc = op(acc,b)
acc = op(acc,c)
acc = op(acc,d)
acc = op(acc,e)

The most important difference between foldLeft and reduceLeft is that the result of foldLeft, and its accumulator, might be of a different type than the elements the operation is applied upon. With reduceLeft, once you applied acc = op(a,b), if acc where of a different type, then op(acc,c) would be a different method. It might work with some implicits or subtypes, but not for entirely different types.

Now, I’m not using the range for nothing, except counting the number of times I want to multiply the matrix. Acc is initialized with “unitMatrixInt(rows)”, which would be a method producing a “unit” matrix of the appropriate size. The “1” of matrices. Then I keep multiplying “this” against this accumulator.

I could create a list of copies of “this”, and then multiply them with reduceLeft too. It would go like this:

def **(n: Int) = List.fill(n)(this) reduceLeft (_*_)

Much simpler, but I wanted to mention foldLeft somehow. :-)

While simplicity has its merits, there is a simple optimization when doing exponentiation by taking advantage of common factors. I went for it in my final implementation:

  def **(n: Int): MatrixInt =
    if (rows != cols  n < 0 =""> unitMatrixInt(rows)
      case 1 => this
      case 2 => this * this
      case odd if odd % 2 == 1 => this ** (odd - 1) * this
      case even => this ** (even / 2) ** 2
    }

I'll make a note on the method name, here. I chose ** as this is one of the common names for such an operation on computer languages. Another, less common, name is ^. I prefer the later, but it has a very low precedence. In fact, this very method would be broken as written, because "(odd - 1) * this" would have precedence over "this ** (odd - 1)", which would completely invert the logic and cause a compilation error.

Now, aside from , ^, &, <>, = and !, :, +, -, * and / and % (ordered by precedence, by the way), all other special characters have a higher precedence. If I called it "@", for instance, then "a * b @ c" would behave as one would expect it too. Unfortunately, no one would expect "@" to be an exponentiation character. But if I were to make extensive use of this library in formulas and expressions, I might well go for something like that, to prevent unexpected errors due to precedence.

Anyway, we need to define unitMatrixInt too. We’ll do so in a companion object. This way, users of this package can create unix matrices if they want to. To create this method, we’ll need a more general factory method, though. We already have most of what we need in createArrays.

But before I get into that, I want to revise my class definition. Presently, this is what we have:

class MatrixInt(array: Array[Array[Int]]) {
  private val self = clone(array)
}

I receive a matrix, and then replicate it internally, so I can be sure it won’t be changed outside my class – and I take care of not changing it inside. But there is a downside to it: “array” continues to be part of my object (see footnote 1)! It becomes a private reference inside my object, and this has some unfortunate consequences. For one thing, since each instance of my class makes reference to the arrays used to build them, those arrays never get deallocated. Now think about all references to “new MatrixInt” inside the class. I use that for every operator. That means I keep TWO copies of a MatrixInt’s array in memory, even when the copy used for construction gets immediately discarded.

This presents a conundrum.On one hand, if the instances of my class are immutable, then I need the full array to begin with. On the other hand, if I receive an array from an external source, I need to copy it before using, as the original array might get changed by whoever passed it to me.

We solve this with the help of a private constructor and an object companion. A private constructor is, of course, a constructor that cannot be accessed from outside the class. The private constructor thing is easy:

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

Now, how to access it? We cannot define another constructor inside MatrixInt with the same signature – receiving an array of array of Int as parameter – because it would conflict with the first. What we do, then, is use an object companion.

Simply put, an object companion is a singleton object – declared with the “object” keyword – which has the same name as a class. That grants both the object and the class visibility of each other’s private definitions.

We’ll then define an “apply” method in that object companion. As you may remember, the “apply” method receives special treatment from the Scala compiler, which enables you to do away with the name – apply – and just pass the parameters to the instance of your object. Since we have to call our object companion “MatrixInt”, then MatrixInt(…) will be converted by Scala into MatrixInt.apply(…). It looks like the instantiation of a new object, only it lacks the “new” keyword.

Let’s define it, then. First, let’s remove the “val self” from the class, and just receive it in its private constructor. As that constructor will only be accessed through MatrixInt’s object companion, that will be safe.

class MatrixInt private (self: Array[Array[Int]]) {
  import MatrixInt._

The import statement makes it possible to access the methods of the object companion without preceding them with “MatrixInt.”. Being a companion object makes the private methods accessible, but does not bring them into scope automatically.

Now, for the object companion, let’s create it, the factory for MatrixInt, and a couple of helpful factories. We’ll also transfer the “clone” and “createArrays” methods to it, as this method does not make reference to a MatrixInt’s own data, so it doesn’t need to be part of each instance. The definition of toArray needs “MatrixInt.” prepended to “clone” to avoid an ambiguous reference, but nothing else should change.

object MatrixInt {
  import MatrixInt._
  implicit def toMatrixInt(m : Array[Array[Int]]) = MatrixInt(m)

  def apply(array: Array[Array[Int]]): MatrixInt = new MatrixInt(clone(array))
  def apply(rows: Int, cols: Int): MatrixInt = MatrixInt(rows, cols, 0)
  def apply(rows: Int, cols: Int, value: Int): MatrixInt =
    new MatrixInt(createArrays(rows, cols, ((_,_) => value)))
  def apply(rows: Int, cols: Int, f: (Int,Int) => Int): MatrixInt =
    new MatrixInt(createArrays(rows, cols, f))
 
  def unitMatrixInt(n: Int) = MatrixInt(n, n, (row,col) => (if (row == col) 1 else 0))

  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))
}

We can now create a MatrixInt by using MatrixInt(Array(Array(…))), but we didn’t – and won’t – change that in the class. The companion object factory produces a clone of the array, just as we wanted. But the arrays produced by MatrixInt’s methods don’t get exposed, so there is no risk in using them. So, by keep using the private constructor inside MatrixInt itself, we prevent unnecessary copies of the arrays produced by each of MatrixInt’s operators.

We’ll finish today with Fibonacci, but we’ll continue later on. I’ll add some methods for completeness to MatrixInt, and then look into MatrixDouble and solve some linear equations with it.

After that, we’ll briefly explore ways to make our Matrix class parametrizable, though I’m still in doubt about performance. Still, I’m reading the comments, and I plan to address them one way or another in these posts.

So, the Fibonacci!

def fibonacci(n: Int) = (MatrixInt(Array(Array(1,1),Array(1,0))) ** n)(1,0)

Footnote 1:

A reader has brought to my attention that, actually, "array" is not being preserved. I haven't found in Scala specification anything about this, so I'm putting it down to optimization. But because it is not explicitly documented, I advise against taking advantage of it.

Do note, please, that if you declare "array" to be protected, val or var, or if any method makes use of it, then it will be preserved in the object.

Matrices 2

I’ll continue today with a little series for beginners. On my first post I forgot something important. One of the very first things I do with a class is create a toString method. This way, I get instant feedback on how my class is doing, as I paste it into the REPL, and try out the methods I defined.

So, let’s think of a toString method. Since all data in my class is stored in the private val “self”, we could simply delegate the task of producing a string:

override def toString = self.toString

That might do at first, but it would be unreasonable to expect something better formatted than that. First, let’s see if we can’t get the output in multiple lines, using mkString:

override def toString = self mkString “\n”

Unfortunately, the output of that looks like this:

res6: String =
[I@196a753
[I@1c36f46

This wouldn’t happen with a List, but an Array is implemented to be compatible with Java, and Java had Arrays have some intrinsic limitations, as they were present in the language from the very beginning, much before generics, and some design decisions about them proved to be unfortunate. Anyway, let’s fix that:

override def toString = self map (_ mkString “ “) mkString “\n”

Better, but we can improve on that still. For one thing, it would be much better if the columns were aligned. Now, to align the columns, we first must know how big a column can get. Or, in other words, we must know the size of the number with the biggest string representation in the matrix.

So, summing it up, we have to look through all of the matrix elements, get their sizes, and choose the biggest one. Getting a number’s size is easy: number.toString.size. Now, if we had a list, this would be easy, but we have an array of arrays. When we have nested collections and we want to unnest them, we use flatMap. For instance:

scala> Array(Array(15,0),Array(2,-30)) flatMap (x => x)
res19: Array[Int] = Array(15, 0, 2, -30)

The flatMap operation we used, (x => x), simply returns what it was given. If all you want is unnest the collections, this suffices. But, actually, we want the number sizes. Since “x” is an Array, not a number, we can’t do flatMap (x => x.toString.size). But we can replace the numbers with their sizes, and then unnest the arrays, like this:

scala> Array(Array(15,0),Array(2,-30)) flatMap (_ map (_.toString.size))
res20: Array[Int] = Array(2, 1, 1, 3)

Next, we want to get the biggest one. We have a collection of number sizes, and we want to reduce it to the single biggest one. We have two such methods in Scala: reduceLeft and reduceRight. Since “max” is commutative, both will produce the same result. But reduceLeft has a better performance, as it iterates over the collection in the order given. The final code should look like this, then:

val maxNumSize = self flatMap (_ map (_.toString.size)) reduceLeft (_ max _)

Pretty easy, eh? My first version had three lines, and was much more confusing to read. It helps stating what your problem is – you’ll often see you can do exactly that in Scala. So, we have a size, what do we do with it? Well, let’s just make sure all numbers are printed to that size. If we format them with “%2d”, for instance, every number will take at least two characters, with the left side being filled with a space as needed. To produce such formatted numbers, we’ll make use of RichString’s “format” method:

override def toString = {
 val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
 val numFormat = "%"+maxNumSize+"d"
 self map (_ map (numFormat format _) mkString “ “) mkString “\n”
}

That’s it! I added a “projection” there to make it a bit more performatic, as it doesn’t produce intermediate arrays. On Scala 2.8, that method gets replaced with “view”, with “projection” remaining as deprecated. Well, that was so easy I want to do more. First, I want to produce little vertical bars on the sides, to make clear the matrix boundaries, in case we are printing them side by side. That just requires changing the inner mkString to mkString(“| “, “ “, “ |”). But I also want to bend the bars up and down, to indicate the top and bottom lines, as well as produce something special for one-line matrices.

For all of them, we are just making small changes in the inner mkString. We have one for the top line, one for the intermediate ones, and one for the bottom line. We could make a list of all these lines, and then use mkString on the list, but we have to be careful, because the “bottom” line and the “intermediate” lines might not exist, and the top line has two different formats depending on whether the matrix has more than one line or not. Taking that into account, here is a much larger, but, hopefully, nicer toString method:

 override def toString = {
   val maxNumSize = self.projection flatMap (_ map (_.toString.size)) reduceLeft (_ max _)
   val numFormat = "%"+maxNumSize+"d"
   def formatNumber(n: Int) = numFormat format n
   val top = rows match {
     case 1 => _row(0) map formatNumber mkString ("< ", " ", " >")
     case _ => _row(0) map formatNumber mkString ("/ ", " ", " \\")
   }
   val middle = rows match {
     case 1 | 2 => Nil
     case _ => self.toList.tail.init map (_ map formatNumber mkString("| "," "," |"))
   }
   val bottom = rows match {
     case 1 => Nil
     case _ => List(_row(rows - 1) map formatNumber mkString ("\\ ", " ", " /"))
   }
  
   top :: middle ::: bottom mkString "\n"
 }

In a parting note, someone suggested one technique to be able to parametrize the class. There are many, but they are slower than using "AnyVal" types directly. If I happen to be mistaken, please do post about it at length!

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)))

Catching Exceptions

A fellow blogger was complaining about Google Code returning an exception instead of null, which made his code look ugly. Not so, says I! Scala (and monads, I suppose) come to rescue:

def exceptionToLeft[T](f: => T): Either[java.lang.Throwable, T] = try {
  Right(f)
} catch {
  case ex => Left(ex)
}


exceptionToLeft(obj.methodWhichMightThrow()) match {
  case Right(x) => // do stuff with x
  case Left(ex) => println("Oops, got exception " + ex.toString)
}

Enjoy!

Update:

Scala 2.8 has library help for this, at scala.util.control.Exception. It does require one to specify the exception or exceptions to be caught, which is a good thing, and it has many more usage modes. You can just ignore the exception for computations returning Unit, you can get the result as Some(T) and transform the exceptions in None, etc.