Algorithm and Data Representation

Dependence of algorithms on the data representation

Traditionally, an algorithm is written using some data structure
I.e.: the algorithm depends on the data structure

Example: the Jacobi algorithm for solving linear equations:

A x = b | write: A = M - N => (M - N) x = b <=> Mx - Nx = b <=> Mx = Nx + b <=> x = M^-1Nx + M^-1b <=> x = Px + Qb

Jacobi Iteration: xⁿ⁺¹ = Pxⁿ + Qb

There are two operations in the Jacobi iteration:

matrix-vector multiplication:
vector-vector addition:

The Jacobi iteration is in principle only dependent on the matrix-vector multiplication and vector-vector addition operations.
However, in practice, when we write the computer program , the algorithm will depend on HOW we represent a matrix and a vector
The algorithm will be different when the matrix and the vector are stored differently

Jacobi Iteration using "Dense" Matrix representation

In the dense representation, the matrix is stored as two-dimensional array

In that case, the matrix multiplcation algorithm is implemented in this way:

class Vector { public: double x[3]; Vector operator+(Vector &v) { Vector z; int i; for (i = 0; i < 3; i = i + 1) z.x[i] = x[i] + v.x[i]; return(z); } }

class DenseMatrix { public: double A[3][3]; Vector operator*(Vector &v) { Vector z; int i, j; for (i = 0; i < 3; i=i+1) z.x[i] = 0; for (i = 0; i < 3; i=i+1) for (j = 0; j < 3; j=j+1) z.x[i] = z.x[i] + A[i][j]*v.x[j]; return z; } };

The Jacobi Iteration algorithm can be written as follows:

double distance(Vector a, Vector b) { double sum = 0; int i; for (i = 0; i < 3; i = i + 1) sum = sum + fabs(a.x[i] - b.x[i]); return(sum); } Vector Jacobi(DenseMatrix & P, DenseMatrix & Q, Vector b) { Vector v1, v2; int i; double d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); }

Notice that

v2 = P*v1 + Q*b;

will invoke the Matrix-vector multiply function defined in the class DenseMatrix

Example Program: (Demo above code)
- Prog file: click here

Jacobi Iteration using a "Sparse" Matrix representation

Recall that Sparse matrices can be stored using the coordinate-wise system... click here

Example:

[ 11 12 14 ] [ 22 23 25 ] [ 31 33 34 ] [ 42 45 46 ] [ 55 ] [ 65 66 67 ] [ 75 77 78 ] [ 87 88 ]

Representation using Coordinate-wise method:

Index 0 1 2 3 4 5 6 7 8 9 10 ----------------------------------------------------------- Val 11 12 14 22 23 25 31 33 34 42 45 Row 0 0 0 1 1 1 2 2 2 3 3 Col 0 1 3 1 2 4 0 2 3 1 4 Index 11 12 13 14 15 16 17 18 19 20 21 ------------------------------------------------------------ Val 46 55 65 66 67 75 77 78 87 88 - Row 3 4 5 5 5 6 6 6 7 7 - Col 5 4 4 5 6 4 6 7 6 7 -

The matrix-vector multiplication using the Coordinate-wise method is implemented as follows:

for (k = 0; k < N; k = k + 1) result[i] = 0; for (k = 0; k < nnz; k = k + 1) result[Row[k]] = result[Row[k]] + Val[k]*d[Col[k]];

When using the Coordinate-wise method , the matrix multiplcation is implemented in this way:

class Vector { public: double x[3]; Vector operator+(Vector &v) { Vector z; int i; for (i = 0; i < 3; i = i + 1) z.x[i] = x[i] + v.x[i]; return(z); } }

class SparseMatrix { public: int nnz; double val[6]; int row[6]; int col[6]; Vector operator*(Vector &v) { Vector z; int k; for (k = 0; k < 3; k = k + 1) z.x[k] = 0; for (k = 0; k < nnz; k = k + 1) z.x[row[k]] = z.x[row[k]] + val[k]*v.x[col[k]]; return z; } };

In order to solve a linear equation using sparse matrices represented by the Coordinate-wise method, we must write another Jacobi function:

Vector Jacobi(SparseMatrix & P, SparseMatrix & Q, Vector b) { Vector v1, v2; int i; double d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); }

Notice that

v2 = P*v1 + Q*b;

will invoke the Matrix-vector multiply function defined in the class SparseMatrix

We CANNOT the first Jacobi function:

Vector Jacobi(DenseMatrix & P, DenseMatrix & Q, Vector b) { ... }

because the input parameters are of different type and therefore, different functions are invoked.

Example Program: (Demo above code)
- Prog file: click here

A more elegant approach...

A more elegant approach is a solution where the algorithm is independent from the underlying representation of the data
This independence is achieved by to adapting the algorithm to the input type

The idealized solution looks as follows:

Vector Jacobi(Matrix P, Matrix Q, Vector b) { Vector v1, v2; int i; double d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); }

This function that implements the Jacobi algorithm using a generic matrix-vector multiplication (and vector-vector addition) operations.

When solving equation in dense representation we write:

DenseMatrix P1, Q1; Vector sol, b; sol = Jacobi( P1, Q1, b);

When solving equation in sparse representation we write:

SparseMatrix P1, Q1; Vector sol, b; sol = Jacobi( P1, Q1, b);

Nothing else need to be changed !!!

This is possible with object oriented programming techniques ...
But before we can explain the technique, we need some background knowledge on derived classes and polymorphism.