Bird

What Gauss-Seidel is Really Doing

I’ve been reading Alan Sokal’s lecture notes “Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms” today. Once I learned to take the word “Hamiltonian” and mentally substitute “function to be minimized”, they are very clearly written.

Anyway, the notes give an explanation of the Gauss-Seidel iterative method for solving linear systems that is so clear, I feel a little cheated that it was never explained to me before. Let’s start with the typical (confusing!) explanation of Gauss-Seidel (as in, say, Wikipedia). You want to solve some linear system

,

for . What you do is decompose into a diagonal matrix , a strictly lower triangular matrix , and a strictly upper triangular matrix , like

.

To solve this, you iterate like so:

.

This works when is symmetric positive definite.

Now, what the hell is going on here? The first observation is that instead of inverting , you can think of minimizing the function

.

(It is easy to see that the global minimum of this is found when .) Now, how to minimize this function? One natural way to approach things would be to just iterate through the coordinates of , minimizing the function over each one. Say we want to minimize with respect to . So, we are going to make the change . Thus, we want to minimize (with respect to d), this thing:

Expanding this, taking the derivative w.r.t. , and solving gives us (where is the unit vector in the th direction)

,

or, equivalently,

.

So, taking the solution at one timestep, , and solving for the solution at the next timestep, we will have, for the first index,

.

Meanwhile, for the second index,

.

And, more generally

.

Now, all is a matter of cranking through the equations.

Which is exactly the iteration we started with.

The fact that this would work when is symmetric positive definite is now obvious– that is when the objective function we started with is bounded below.

The understanding that Gauss-Seidel is just a convenient way to implement coordinate descent also helps to get intuition for the numerical properties of Gauss-Seidel (e.g. quickly removing "high frequency error", and very slowly removing "low frequency error").

Numbers in bar graphs

I spent far too much time today figuring out how to put numbers in bar graphs. i.e. how to transform this:

into this:

When reading papers, I often wish numbers were included like this, and so resolved a while ago to always do so myself.  Plotting the numbers allows someone else to replot them when doing some sort of a comparison.  (With out needing to implement your algorithm, or ask you for your doubtlessly long-lost data, or (ugh) physically measuring the lengths of your bars.)  However, now seeing how tedious this can be, I understand why it is rarely done.

In any case, I wrote a function add_bar_nums.m that you can run after doing a bar plot that will add the numbers like above.  There are options for rotating the text, changing the display format, etc.

A simple explanation of reverse-mode automatic differentiation

My previous rant about automatic differentiation generated several requests for an explanation of how it works. This can be confusing because there are different types of automatic differentiation (forward-mode, reverse-mode, hybrids.) This is my attempt to explain the basic idea of reverse-mode autodiff as simply as possible.

Reverse-mode automatic differentiation is most attractive when you have a function that takes inputs , and produces a single output . We want the derivatives of that function, , for all .

Point #1: Any differentiable algorithm can be translated into a sequence of assignments of basic operations.

Forward-Prop

for

Here, each function is some very basic operation (e.g. addition, multiplication, a logarithm) and denotes the set of “parents” of . So, for example, if and , then .

It would be extremely tedious, of course, to actually write an algorithm in this “expression graph” form. So, autodiff tools create this representation automatically from high-level source code.

Point #2: Given an algorithm in the previous format, it is easy to compute its derivatives.

The essential point here is just the application of the chain rule.

Applying this, we can compute all the derivatives in reverse order.

Back-Prop

for

That’s it!  Just create an expression graph representation of the algorithm and differentiate each basic operation in reverse order using calc 101 rules.

Other stuff:

• No, this is not the same thing as symbolic differentiation.  This should be obvious:  Most algorithms don’t even have simple symbolic representations.  And, even if yours does,  it is possible that it “explodes” upon symbolic differentiation.  As a contrived example, try computing the derivative of .

• The complexity of the back-prop step is the same as the forward propagation step.
• In machine learning, functions from N inputs to one output come up all the time:  The N inputs are parameters defining a model, and the 1 output is a loss, measuring how well the model fits training data.  The gradient can be fed into an optimization routine to fit the model to data.
• There are two common ways of implementing this:

1. Operator Overloading.  One can create a new variable type that has all the common operations of numeric types, which automatically creates an expression graph when the program is run.   One can then call the back-prop routine on this expression graph.  Hence,  one does not need to modify the program, just replace each numeric type with this new type.  This is fairly easy to implement, and very easy to use.  David Gay’s RAD toolbox for C++ is a good example, which I use all the time.
   The major downside of operator overloading is efficiency:  current compilers will not optimize the backprop code.  Essentially, this  step is interpreted.  Thus, one finds in practice a non-negligible overhead of, say, 2-15 times the complexity of the original algorithm using a native numeric type.  (The overhead depends on how much the original code benefits from compiler optimizations.)
2. Source code transformation. Alternatively, one could write a program that examines the source code of the original program, and transforms this into source code computing the derivatives.  This is much harder to implement, unless one is using a language like Lisp with very uniform syntax.  However, because the backprop source code produced is then optimized like normal code, it offers the potential of zero overhead compared with manually computed derivatives.

• If it isn’t convenient to use automatic differentiation, one can also use “manual automatic differentation”.  That is, to compute the derivatives, just attack each intermediate value your algorithm computes, in reverse order.
• Some of the most interesting work on autodiff comes from Pearlmutter and Siskind, who have produced a system called Stalingrad for a subset of scheme that allows for crazy things like taking derivatives of code that itself is taking derivates.  (So you can, for example, produce Hessians.)  I think they wouldn’t mind hearing from potential users.

Prediction markets, monte carlo, and loss functions

There has been some discussion lately about how to evaluate the performance of different prediction markets (like Intrade), and predictors (like Nate Silver) at guessing the winners of elections, or Oscars. Who is making the best predictions? If everyone simply made a guess for the winner of each state or award, we could evaluate performance easily: whoever guesses the most outcomes correctly is making the best predictions. But what do we do if the predictors provide us full probabilities of the different outcomes? Intuitively, someone who gives 99% probability of an event that doesn’t occur is much “more wrong” than someone who gives only a 51% probability.

                        Nate    Intrade     Winner
Best Picture
Slumdog                .990      .903        X
Milk                   .010      .040
Frost/Nixon                      .013
Benjamin Button                  .080
The Reader                       .030

Best Director
Danny Boyle            .997      .900        X
Gus Van Sant           .001      .059
David Fincher          .001      .050
Ron Howard
Steven Daldry           

Lead Actress
Kate Winslet           .676      .850        X
Meryl Streep           .324      .150
Anne Hathaway                    .044
Melissa Leo
Angelina Jolie         

Lead Actor
Mickey Rourke          .711      .700
Sean Penn              .190      .335        X
Brad Pitt              .059
Frank Langella         .034      .049
Richard Jenkins        .005

Supporting Actress
Taraji P. Henson       .510      .190
Penélope Cruz          .246      .588        X
Viola Davis            .116      .199
Amy Adams              .116
Marisa Tomei           .012

Supporting Actor
Heath Ledger           .858      .950        X
Josh Brolin            .050      .050
Philip Seymour Hoffman .044
Michael Shannon        .036
Robert Downey Jr.      .012

Let us think about this situation from the perspective of “bent coin predictors”. Let’s say we have a pool of 100 bent coins, each of which has some unknown probability of ending up heads. We have a number of people who reckon they can estimate that probability by looking at the coin. Denote the prediction of guesser for coin by After predictions have been made, we flip all the coins. Now, how do we find the best guesser?

What we want, is to measure how close is to . Since we don’t know , this seems impossible. In some sense, it is impossible, but let us fantasize for a moment. Suppose that instead of all the coin flips, someone actually revealed all the true probabilities to us. Then, what would we do? There is no single best answer. One reasonable way to measure the quality of the guess would be the sum of squares difference

or, equivalently,

.

MONTE-CARLO APPROXIMATIONS

Now, of course, we can’t calculate either of the above quantities. We only have a single result from flipping each coin. The central idea here is that we can use what is known as a Monte-Carlo approximation. This is a very simple idea. Suppose we would like to calculate the expected value of some function .

Now, suppose that we don’t know , but we can simulate . That is, by running some sort of experiment, we can get some random value , whose probability is . If we draw many such values, then we can approximate the above by:

It would be interesting to look at how accurate various predictors were for the 2008 election from this perspective.

As an example, suppose we want to know the average amount that a slot machine pays out. We could approximate this by playing the machine 1000 times, and calculating the average observed payout.

LOSS FUNCTIONS

How can we apply loss functions to prediction markets? Notice that we can make the following simplification

.

The first term, is hard to estimate, but fortunately we don’t need to, because it doesn’t depend on the predictions. If we ignore this term for all guessers, we still have a valid relative rating of the guessers.

As for the second term, we can apply the Monte-Carlo technique from above. Let denote the outcome of flipping coin . We make the approximation

.

This is a very noisy approximation. However, it is unbiased: If we average over many coins, we will get something very close to the true value. (This works exactly the same way that we could approximate the average payout of slot machines in a casino by playing 1000 random machines and then averaging the payouts.)

Happily, the final term,

,

can be computed exactly, since the guessers have provided .

THE OSCARS

Now, let’s apply the above theory to the Oscar predictions.

Taking a zero for the empty entries, and normalizing each prediction, we obtain the scores:

Quadratic loss:

538:     -0.6235
intrade: -0.7925

Remember, less is better, so this is a clear win for intrade.

OTHER LOSS FUNCTIONS

Another reasonable way to measure errors would be the KL-divergence between and

.

Again, dropping a constant term, we get the loss

.

For historical reasons, let’s call this this “conditional likelihood” loss.

Conditional likelihood loss:

538:     0.6032
intrade: 0.3699

Again, intrade looks much better. Then again, of course, maybe intrade was just lucky. Five predictions isn’t a large number to average over.

I’d love to see this type of analysis applied to the state-by-state results of the 2008 elections.

Matlab code (probably also works with Octave) is here.

The Stalin Compiler

Stalin is a questionably named Scheme compiler written by Jeffrey Mark Siskind that can supposedly create binaries as fast or faster than Fortran or C for numerical problems. To test this, I tried creating a simple program to numerically integrate from 0 to 10000. To make things interesting, I used a manual Newton’s method implementation of sqrt from SICP.  The integration is done by a simple tail-recursive method.

The scheme code is very pretty:

(define (sqrt-iter guess x)
  (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

(define (improve guess x)
  (average guess (/ x guess)))

(define (average x y)
  (/ (+ x y) 2))

(define (good-enough? guess x)
  (< (abs (- (* guess guess) x)) 0.001))

(define (mysqrt x)
  (sqrt-iter 1.0 x))

(define (int x acc step)
  (if (>= x 10000.0)
      acc
      (int (+ x step)
           (+ acc (* step (mysqrt x)))
           step)))

(write (int 0.0 0.0 .001))

I then converted this to C. It is pretty much a transliteration, except it uses a for loop instead of recursion to accumulate the values:

#include "stdio.h"

double improve(double guess, double x);
double average(double x, double y);

double sqrt_iter(double guess, double x){
  if( good_enough(guess, x))
    return guess;
  else
    return sqrt_iter( improve(guess,x), x);
}

double improve(double guess, double x){
  return average(guess, x/guess);
}

double average(double x, double y){
  return (x+y)/2;
}

int good_enough(double guess, double x){
  if (abs(guess*guess-x)<.001)
    return 1;
  return 0;
}

double mysqrt(double x){
  return sqrt_iter(1.0, x);
}

main(){
  double rez = 0;
  double x;
  double step = .001;
  for(x=0; x<= 10000; x+=step)
    rez += mysqrt(x)*step;
  printf("%f\n", rez);
}

I compiled the two methods via:

stalin -On -d -copt -O3 int.sc
gcc -O3 int.c

The results are:
Stalin: 1.90s
gcc: 3.61s

If you declare every method inline in C, you get:

gcc-inline: 3.28s

For reference, I also compiled the code with chicken scheme, using the -optimize-level 3 compiler flag.

Chicken Scheme: 27.9s.

Some issues:

• The Scheme code is more appealing on the whole, though it suffers from a lack of infix notation: (< (abs (- (* guess guess) x)) 0.001) is definitely harder to read than abs(guess*guess-x)<.001.  I wonder if more familiarity with prefix code would reduce this.
• All three methods give slightly different results, which is worrisome.  In particular, Stalin is correct to 6 digits, whilst gcc and chicken are correct to 7.
• Stalin apparently does not include macros.  One would think it would be easy to use the macro system from a different scheme compiler, and then send the source to stalin for final compilation.
• Stalin is extremely slow to compile.  In principle this isn’t a big deal: you can debug using a different scheme compiler.  Still, Stalin seems to be somewhat less robust to edge cases, than at least chicken scheme.
• It is amazing that Scheme code with no type declarations can beat C by almost a factor of 2.
• Though in principle Stalin produces intermediate c code, it is utterly alien and low-level.  I have not been able to determine exactly what options Stalin is using when it calls gcc on the source code.  That could account for some of the difference.
• A detailed writeup on Stalin is here.

Automatic Differentiation: The most criminally underused tool in the potential machine learning toolbox?

I recently got back reviews of a paper in which I used automatic differentiation.  Therein, a reviewer clearly thought I was using finite difference, or “numerical” differentiation. This has led me to wondering: Why don’t machine learning people use automatic differentiation more?  Why don’t they use it…constantly? Before recklessly speculating on the answer, let me briefly review what automatic differentiation (henceforth “autodiff”) is. Specifically, I will be talking about reverse-mode autodiff.

(Here, I will use “subroutine” to mean a function in a computer programming language, and “function” to mean a mathematical function.)

It works like this:

1. You write a subroutine to compute a function .  (e.g. in C++ or Fortran).  You know to be differentiable, but don’t feel like writing a subroutine to compute .
2. You point some autodiff software at your subroutine.  It produces a subroutine to compute the gradient.
3. That new subroutine has the same complexity as the original function!
   1. It does not depend on the dimensionality of .
4. It also does not suffer from round-off errors!

To take a specific example, suppose that you have written a subroutine to evaluate a neural network.  If you perform autodiff on that subroutine, it will produce code that is equivalent to the backpropagation algorithm. (Incidentally, autodiff significantly predates the invention of backprop).

Why should autodiff be possible?  It is embarasingly obvious in retrospect:  Complex subroutines consist of many elementary operations.  Each of those is differentiable.  Apply the calc 101 chain rule to the expression graph of all these operations.

A lot of papers in machine learning (including many papers I like) go like this:

1. Invent a new type of model or a new loss function
2. Manually crank out the derivatives.  (The main technical difficulty of the paper.)
3. Learn by plugging the derivatives into an optimization procedure.  (Usually L-BFGS or stochastic gradient.)
4. Experimental results.

It is bizarre that the main technical contribution of so many papers seems to be something that computers can do for us automatically.  We would be better off just considering autodiff part of the optimization procedure, and directly plugging in the objective function from step 1.  In my opinion, this is actually harmful to the field.  Before discussing why, let’s consider why autodiff is so little used:

1. People don’t know about it.
2. People know about it, but don’t use it because they want to pad their papers with technically formidable derivations of gradients.
3. People know about it, but don’t use it for some valid reasons I’m not aware of.

I can’t comment on (3) by definition.  I think the answer is (1), though I’m not sure.  Part of the problem is that “automatic differentiation” sounds like something you know, even if you actually have no idea what it is.  I sometimes get funny looks from people when I claim that both of the following are true.

.

Further evidence for (1) is that many papers, after deriving their gradient, proclaim that “this algorithm for computing the gradient is the same order of complexity as evaluating the objective function,” seeming to imply it is fortunate that a fast gradient algorithm exists.

---

Now, why bother complaining about this?  Are those manual derivatives actually hurting anything?  A minor reason is that it is distracting.  Important ideas get obscured by the details, and valuable paper space (often a lot of space) is taken up for the gradient derivation rather than more productive uses.  Similarly, valuable researcher time is wasted.

In my view a more significant downside is that the habit of manually deriving gradients restricts the community to only using computational structures we are capable of manually deriving gradients for.  This is a huge restriction on the kinds of computational machinery that could be used, and the kinds of problems that could be addressed.

As a simple example, consider learning parameters for some type of image filtering algorithm.  Imaging problems suffer from boundary problems.  It is not possible to compute a filter response at the edge of an image.  A common trick is to “flip” the image over the boundary to provide the nonexistent measurements.  This poses no difficulty at all to autodiff, but borders on impossible to account for analytically.  Hence, I suspect, people simply refrain from researching these types of algorithms.  That is a real cost.

For C++ users, it is probably easiest to get started with David Gay’s RAD toolbox.  (See also the paper.)

Caveats:

• Often, derivatives are needed for analysis, not just to plug into an optimization routine.  Of course, these derivatives should always remain.
• When derivatives are reasonably simple, it can be informative to see the equations.  Oftentimes, however, an algorithm is derived for computing the gradient.  In this cases, it would almost always be easier for someone trying to implement the method to install an autodiff tool than try to implement the gradient algorithm.

Update: Answers to some questions:

Q) What’s the difference between autodiff and symbolic diff?

A) They are totally different. The biggest difference is that autodiff can differentiate algorithms, not just expressions. Consider the following code:

function f(x)
  y = x;
  for i=1...100
    y = sin(x+y);
  return y

Automatic differentiation can differentiate that, easily, in the same time as the original code. Symbolic differentiation would lead to a huge expression that would take much more time to compute.

Q) What about non-differentiable functions?

A) No problem, as long as the function is differentiable at the place you try to compute the gradient.

Q) Why don’t you care about convexity?

A) I do care about convexity, of course.

Hessian-Vector products

You have some function .  You have figured out how to compute it’s gradient, .  Now, however, you find that you are implementing some algorithm (like, say, Stochastic Meta Descent), and you need to compute the product of the Hessian with certain vectors.  You become very upset, because either A) you don’t feel like deriving the Hessian (probable), or B) the Hessian has elements, where is the length of and that is too big to deal with (more probable).  What to do?  Behold:

Consider the Hessian-Vector product we want to compute, .  For small ,

And so,

This trick above has been around apparently forever.  The approximation becomes exact in the limit .  Of course, for small , numerical problems will also start to kill you.  Pearlmutter’s algorithm is a way to compute with the same complexity, with out suffering rounding errors.  Unfortunately, Pearlmutter’s algorithm is kind of complex, while the above is absolutely trivial.

Update: Perlmutter himself comments below that if we want to use the finite difference trick, we would do better to use the approximation:

This expression will be closer to the true value for larger , meaning that we are less likely to get hurt by rounding error. This is nicely illustrated on page 6 of these notes.

Why does regularization work?

When fitting statistical models, we usually need to “regularize” the model. The simplest example is probably linear regression. Take some training data, . Given a vector of weights , the total squared distance is

So to fit the model, we might find to minimize the above loss. Commonly, (particularly when has many dimensions), we find t\hat the above procedure leads to overfitting: very large weights t\hat fit the training data very well, but poorly predict future data. “Regularization” means modifying the optimization problem to prefer small weights. Consider

Here, trades off how much we care about the model fitting well versus how we care about being small. Of course, this is a much more general concept that linear regression, and one could pick many alternatives to the squared norm to measure how simple a model is. Let’s abstract a bit, and just write some function to represent how well the model fits the data, and some function to represent how simple it is.

The obvious question is: how do we pick ? This, of course, has been the subject of much theoretical study. What is commonly done in practice is this:

1. Divide the data into a training set and test set.
2. Fit the model for a range of different s using only the training set.
3. For each of the models fit in step 2, check how well the resulting weights fit the test data.
4. Output the weights that perform best on test data.

(One can also retrain on all the data using the that did best in step 2.)

Now, there are many ways to measure simplicity. (In the example above, we might consider , , , , …). Which one to use? If you drank too much coffee one afternoon, you might decide to include a bunch of regulizers simultaneously, each with a corresponding regularization parameter, ending up with a problem like

.

And yet, staring at that equation, things take a sinister hue: if we include too many regularization parameters, won’t we eventually overfit the regularization parameters? And furthermore, won’t it be very hard to test all possible combinations of to some reasonable resolution? What should we do?

And now I will finally get to the point.

Recently, I’ve been looking at regularization in a different way, which seems very obvious in retrospect. The idea is that when searching for the optimal regularization parameters, we are fitting a model– a model that just happens to be defined in an odd way.

First, let me define some notation. Let denote the training set, and let denote the test set. Now, define

Where is the function measuring how well the model fits the training data. So, for a given regularization parameter, returns the weights solving the optimization problem. Given that, the optimal \lambda will the the one minimizing this equation:

This extends naturally to the setting where we have a vector of regularization parameters .

From this perspective, doing an optimization over regularization parameters makes sense: it is just fitting the parameters to the data – it just so happens that the objective function is implicitly defined by a minimization over the data .

Regularization works because it is fitting a single parameter to some data. (In general, one is safe fitting a single parameter to any reasonably sized dataset although there are exceptions.) Fitting multiple regularization parameters should work, supposing the test data is big enough to support them. (There are computational issues I’m not discussing here with doing this, stemming from the fact that isn’t known in closed form, but is implicitly defined by a minimization. So far as I know, this is an open problem, although I suspect the solution is “implicit differentiation”)

Even regularizing the regularization parameters isn’t necessarily a crazy idea, though clearly isn’t the way to do it. (Simpler models correspond to larger regularization constants. Maybe ?.) Then, you can regularize those parameters!

Well. In retrospect that wasn’t quite so simple as it seemed.

It’s a strong bet that this is a standard interpretation in some communities.

Update: The idea of fitting multiple parameters using implicit differentiation is published in the 1996 paper Adaptive regularization in neural network modeling by Jan Larsen, Claus Svarer, Lars Nonboe Andersen and Lars Kai Hansen.

Mandelbrot in Scala

With my growing frustration with Matlab, I’ve been looking for a while for a language that was

1. Garbage collected
2. Good notation for numerical computation
3. Fast enough for numerical computation

After a long search, I think I’ve finally found my home in Scala. Today I did a very trivial, self-contained computation of some images of the Mandelbrot set to test it out.

When learning Scala, I couldn’t find enough examples of numerical stuff, so I thought I would post this here for those potentially interested. Here’s a first attempt:

import java.io._
import java.lang.Math

object mandelbrot {
  // tiny complex number class including syntactic sugar for basic operations
  class Complex(val a: Double, val b: Double){
    // represents the complex number a + b*i
    def +(that: Complex) = new Complex(this.a+that.a,this.b+that.b)
    def *(that: Complex) = new Complex(this.a*that.a-this.b*that.b,this.a*that.b+that.a*this.b)
    def abs() = Math.sqrt(this.a*this.a + this.b*this.b)
  }

  def run(n: Int, level: Int) : Unit = {
    val out = new FileOutputStream("scalaimage.pgm")
    out.write(("P5\n"+n+" "+n+"\n255\n").getBytes())

    for {y0 <- 0 until n
	 x0 <- 0 until n }{
	   // y0 and x0 are in pixel integers
	   // x  and y  are real number coordinates
	   val x = -2.0 + x0*3.0/n
	   val y = -1.5 + y0*3.0/n

	   var z = new Complex(0,0)
	   var c = new Complex(x,y)
	   for(i <- 0 until level)
	     z = z*z + c 

	   if (z.abs < 2)
	     out.write(0);
	   else
	     out.write(255);
	 }
    out.close()
  }

  def main(args: Array[String]) {
    run(Integer.parseInt(args(1)), Integer.parseInt(args(0)))
  }
}

Not bad, right? Notice how painlessly we create the Complex class, along with the natural syntax for manipulating the numbers. What’s more, the result is very fast. For a 2048 by 2048 image, I get results in less than 30 seconds. (Just consider running this algorithm in Matlab. Ha!)

Save this to a file “mandelbrot1.scala”, and compile and run with:

scalac mandelbrot1.scala
scala mandelbrot 100 2048

The first argument is many times to try the mandelbrot iteration, and the second argument is how big of an image to output. The result is:

(I manually converted the .pgm to .png to save space– click to get the full resolution version.)

A slightly more complex version follows. I added notation for the operators *= and +=, which speeds things up by about 10%. I also now color the pixels by how many iterations it takes the mandelbrot iterations to diverge.

import java.io._
import java.lang.Math

object mandelbrot {
  // tiny complex number class including syntactic sugar for basic operations
  class Complex(var a: Double, var b: Double){
    // represents the complex number a + b*i
    def +(that: Complex) = new Complex(this.a+that.a,this.b+that.b)
    def *(that: Complex) = new Complex(this.a*that.a-this.b*that.b,this.a*that.b+that.a*this.b)
    def abs() = Math.sqrt(this.a*this.a + this.b*this.b)
    def *=(that: Complex) ={
      val newa = this.a*that.a-this.b*that.b
      this.b = this.a*that.b+that.a*this.b
      this.a = newa
      this
    }
    def +=(that: Complex)={
      this.a += that.a
      this.b += that.b
      this
    }
  }

  def run(n: Int, level: Int) : Unit = {
    val out = new FileOutputStream("scalaimage.pgm")
    out.write(("P5\n"+n+" "+n+"\n255\n").getBytes())

    for {y0 <- 0 until n
	 x0 <- 0 until n }{
	   // y0 and x0 are in pixel integers
	   // x  and y  are real number coordinates
	   val x = -2.0 + x0*3.0/n
	   val y = -1.5 + y0*3.0/n

	   var z = new Complex(0,0)
	   var c = new Complex(x,y)
	   var i = 0
	   do {
	     z *= z
	     z += c
	     i += 1
	   } while( z.abs < 2 && i < level)

	   if (z.abs < 2)
	     out.write(0);
	   else
	     out.write( (i*255.0/level).toInt );
	 }
    out.close()
  }

  def main(args: Array[String]){
    run(Integer.parseInt(args(1)), Integer.parseInt(args(0)))
  }
}

Compile and run as before. The output is:

I’m not sure how many people are using Scala for numerical applications, but it looks very good so far. There are a few minor tricks that turn out to be incredibly convenient. For example, if you define any function func for some class a, in a normal language you would call it by

a.func(b).

Scala defines this, but also gives you automatically the notation

a func b.

Adding in Scala’s non-mandatory type declarations (as above, notice that x0 and x aren’t declared to be Int or Double) tricks for implicit conversion between types, and it looks like it is possible to recreate almost all of Matlab’s synactic sugar for matrix manipulation, while actually allowing reasonable speed for non-vectorizable code, and reasonable facilities for abstraction. I hope to soon be an ex-Matlab user.

The scalala library looks very promising, though I haven’t taken the time to install it yet, and it looks like absolutely no one is using it. (dramage, give us some documentation!)