Can't get Jacobi algorithm to work in Objective-C

Question

Hi,

For some reason, I can't get this program to work. I've had other CS majors look at it and they can't figure it out either.

This program performs the Jacobi algorithm (you can see step-by-step instructions and a MATLAB implementation here). BTW, it's different from the Wikipedia article of the same name.

Since NSArray is one-dimensional, I added a method that makes it act like a two-dimensional C array. After running the Jacobi algorithm many times, the diagonal entries in the NSArray (i[0][0], i[1][1], etc.) are supposed to get bigger and the others approach 0. For some reason though, they all increase exponentially. For instance, i[2][4] should equal 0.0000009, not 9999999, while i[2][2] should be big.

Thanks,

Chris

NSArray+Matrix.m

@implementation NSArray (Matrix)

@dynamic offValue, transposed;

- (double)offValue {
    double sum = 0.0;
    for ( MatrixItem *item in self )
        if ( item.nonDiagonal )
            sum += pow( item.value, 2.0 );
    return sum;
}

- (NSMutableArray *)transposed {
    NSMutableArray *transpose = [[[NSMutableArray alloc] init] autorelease];
    int i, j;

    for ( i = 0; i < 5; i++ ) {
        for ( j = 0; j < 5; j++ ) {
            [transpose addObject:[self objectAtRow:j andColumn:i]];
        }
    }
    return transpose;
}

- (id)objectAtRow:(NSUInteger)row andColumn:(NSUInteger)column {
    NSUInteger index = 5 * row + column;
    return [self objectAtIndex:index];
}

- (NSMutableArray *)multiplyWithMatrix:(NSArray *)array {
    NSMutableArray *result = [[NSMutableArray alloc] init];
    int i = 0, j = 0, k = 0;
    double value;

    for ( i = 0; i < 5; i++ ) {
        for ( j = 0; j < 5; j++ ) {
            value = 0.0; // (JeremyP's answer)
            for ( k = 0; k < 5; k++ ) {
                MatrixItem *firstItem = [self objectAtRow:i andColumn:k];
                MatrixItem *secondItem = [array objectAtRow:k andColumn:j];
                value += firstItem.value * secondItem.value;
            }
            MatrixItem *item = [[MatrixItem alloc] initWithValue:value];
            item.row = i;
            item.column = j;
            [result addObject:item];
        }
    }
    return result;
}
@end

Jacobi_AlgorithmAppDelegate.m

// ...

- (void)jacobiAlgorithmWithEntry:(MatrixItem *)entry {
    MatrixItem *b11 = [matrix objectAtRow:entry.row andColumn:entry.row];
    MatrixItem *b22 = [matrix objectAtRow:entry.column andColumn:entry.column];

    double muPlus = ( b22.value + b11.value ) / 2.0;
    muPlus += sqrt( pow((b22.value - b11.value), 2.0) + 4.0 * pow(entry.value, 2.0) );

    Vector *u1 = [[[Vector alloc] initWithX:(-1.0 * entry.value) andY:(b11.value - muPlus)] autorelease];
    [u1 normalize];
    Vector *u2 = [[[Vector alloc] initWithX:-u1.y andY:u1.x] autorelease];

    NSMutableArray *g = [[[NSMutableArray alloc] init] autorelease];
    for ( int i = 0; i <= 24; i++ ) {
        MatrixItem *item = [[[MatrixItem alloc] init] autorelease];
        if ( i == 6*entry.row )
            item.value = u1.x;
        else if ( i == 6*entry.column )
            item.value = u2.y;
        else if ( i == ( 5*entry.row + entry.column ) || i == ( 5*entry.column + entry.row ) )
            item.value = u1.y;
        else if ( i % 6 == 0 )
            item.value = 1.0;
        else
            item.value = 0.0;

        [g addObject:item];
    }

    NSMutableArray *firstResult = [[g.transposed multiplyWithMatrix:matrix] autorelease];

    matrix = [firstResult multiplyWithMatrix:g];
}

// ...

Answer 1

Have you got any unit tests for your matrix category? I mean, are you certain that the multiplication algorithm works? I would say that initialising value to 0 happens in the wrong loop. I think you need to do it inside the j loop.

A couple of other observations:

• You don't need the @dynamic property declaration because you are defining the implementation of the properties yourself.
• Consider creating your own Matrix class that wraps a normal C array of doubles. You might find the implementation a bit simpler.

Answer 2

When you add the square root term to muPlus, you don't divide by two. The calculation should be either:

double muPlus = ( b22.value + b11.value ) / 2.0;
muPlus += sqrt( pow((b22.value - b11.value), 2.0) 
                + 4.0 * pow(entry.value, 2.0) 
              ) / 2.0;

or:

double muPlus = ( b22.value + b11.value );
muPlus += sqrt( pow((b22.value - b11.value), 2.0) 
                + 4.0 * pow(entry.value, 2.0) );
muPlus /= 2.0;

Also, you assign u1.y to both G_r,c and G_c,r. From the algorithm description, you want G_r,c=U_1,2 (or u1.y) and G_c,r=U_2,1 (or u2.x). Note that you don't actually need u2; you can substitute -u1.y for u2.x and u1.x for u2.y.

Off-Topic

According to the Fundamental Rule of Cocoa Memory Management, -[NSArray multiplyWithMatrix:] should return an autoreleased array, since the multiplicand should relinquish ownership. Also, you should use accessors to assign GT * A * G to matrix rather than doing it directly so that it can be properly managed.

Since most of the tests in the loop to fill out g will be false during each iteration, it's most likely more efficient to fill g with some default values and then update g. You could create a zero matrix, then set the diagonal to ones, then fill in the values from U, or you could create an identity matrix (leave the i%6 == 0 test in the loop) then fill in the values from U. Profile each of the three approaches.