Caching policies and techniques for matrices

Question

Hi all, as explained before, I'm currently working on a small linear algebra library to use in a personal project. Matrices are implemented as C++ vectors and element assignment ( a(i,j) = v; ) is delegated to the assignment to the vector's elements. For my project I'll need to solve tons of square equation systems and, in order to do that, I implemented the LU factorization (Gaussian Elimination) for square matrices. In the current implementation I'm avoiding to recalculate each time the LU factorization by caching the L and U matrices, the problem is that since I'm delegating the element assignment to vector, I can't find a way to say if the matrix is being changed and whether to recalculate the factorization. Any ideas on how to solve this?

Thank you

Answer 1

Your description makes it sound like you're doing the LU decomposition in-place. That's certainly more efficient from a memory point of view. That means overwriting the matrix values as you perform the decomposition. If that's true, then "being changed and whether to recalculate the factorization" is a moot point. You lose the original matrix when you overwrite it with the LU decomposition.

In the event that you are NOT overwriting the original, it also sounds like you'd want to recalculate the decomposition whenever you give a matrix element a new value. I'd recommend that you not do that. It seems inefficient to me. If a client wanted to alter many values, they probably wouldn't want to pay the cost of another LU decomposition until they were all done.

You can try a factory interface for matrix transformations/decompositions. It's a simple one that will take in a Matrix and return a (decomposed) matrix. You get to keep your original matrix that way; the return value is a new instance. You can change the original and then pass it to the factory to recalculate the LU decomposition. It costs you memory, which can be a problem for very large matricies.

In Java, I'd write it like this:

public interface MatrixDecomposition
{  
    Matrix decompose(Matrix original);
}

In C++, it'd be a pure virtual function. (Been too long - I can't remember the syntax.)

There are other types of decomposition (e.g., QR, SVD, etc.), so this design will nicely accomodate those when you need them. Just write another implementation for the interface and Bob's your uncle.

Lots of physics problems are characterized by "sparse" matricies, which have a bandwidth of non-zero values clustered around the diagonal and zeroes outside. If you use special techniques that don't store the zero values outside the bandwidth you can solve larger problems in memory.

Answer 2

If I understand correctly you need to check during the execution of your code whether a matrix has changed or not.

Well, vectors don't support such functionality. However, what you can do is write a Matrix class of your own, add such functionality to it and use it instead of vectors.

An example implementation could be:

class Matrix {
public:
    Matrix() : hasChanged(false) {}

    double setElement(int i, int j, double value) {
        innerStorage[i][j] = value;
        hasChanged = true;
    }

    double getElement(int i, int j) {
        return innerStorage[i][j];
    }

    void clearHasChangedFlag() {
        hasChanged = false;
    }

private:
    vector<vector<double> > innerStorage;
    bool hasChanged;
}

Answer 3

The setElement and getElement version suggested by Frederick was an option, but it's just too boring to use it. I want to access a(i,j) in reading and writing with the same syntax, withouth bothering if I'm modifying the matrix or not. Probably the best thing to do is let the user modify the Matrix and delegate to him the responsibility to recalculate the LU factorization when he judges that's necessary. Yet it will be very boring.

Answer 4

What about keeping a complete hidden copy of the Matrix used with the last LU factorization and check it in O(n*m) time before re-doing it? Ugly but could work.

Answer 5

template<class T>
class matrix {
public:
 class accessor {
 public:
  accessor(T& dest, matrix& parent) : dest(dest), parent(parent) { }
  operator T const& () const { return dest; }
  accessor& operator=(T const& t) { dest = t; parent.invalidate_cache(); return *this; }
 private:
  T& dest;
  matrix& parent;
 };

 // replace those with actual implementation..
 accessor operator()(int x, int y) {
  static T t; return accessor(t, *this);
 }
 T const& operator()(int x, int y) const {
  static T t; return t;
 }

private:
 void invalidate_cache() { cout << "Cache invalidated !!\n"; }
 vector<T> impl;
};

thanks go to to ##iso-c++ @ irc.freenode.net for some helpful corrections