Generic matrix library in Java

Question

I have this repository. The most important source files follow.

com.github.coderodde.math.linear.matrix.AbstractMatrix.java:

package com.github.coderodde.math.linear.matrix;

import java.util.Objects;

/**
 * This abstract class defines the API for the matrix data types.
 * 
 * @param <M> the actual implementing matrix type.
 * @param <E> the matrix element type.
 * 
 * @author Rodion "rodde" Efremov
 * @version 1.6 (Aug 13, 2023)
 * @since 1.6 (Aug 13, 2023)
 */
public abstract class AbstractMatrix<M extends AbstractMatrix<M, E>, E> {
    
    /**
     * The width of this matrix.
     */
    protected final int width;
    
    /**
     * The height of this matrix.
     */
    protected final int height;
    
    /**
     * The field element API.
     */
    protected final FieldElement<E> fieldElements;
    
    protected AbstractMatrix(int width, 
                             int height, 
                             FieldElement<E> fieldElements) {
        
        this.width = checkWidth(width);
        this.height = checkHeight(height);
        this.fieldElements = Objects.requireNonNull(fieldElements);;
    }
    
    public int getWidth() {
        return width;
    }
    
    public int getHeight() {
        return height;
    }
    
    /**
     * Sets each element in this matrix to its negative.
     */
    public abstract void negate();
    
    /**
     * Returns a copy of this matrix with each element negated. After this 
     * operation, this matrix remains intact.
     * 
     * @return new, negated matrix.
     */
    public abstract M immutableNegate();
    
    /**
     * Adds the input matrix to this matrix.
     * 
     * @param other the matrix to add to this matrix. 
     */
    public abstract void add(M other);
    
    /**
     * Returns a copy of this matrix with elements from {@code other} added to 
     * it. After this operation, this matrix remains intact.
     * 
     * @param other the matrix to add.
     * 
     * @return copy of this matrix with input elements added.
     */
    public abstract M immutableAdd(M other);
    
    /**
     * Creates a new matrix and sets it to the product of this and {@code right}
     * matrices. After this operation, this matrix remains intact.
     * 
     * @param right the right hand matrix in the product. This matrix is the 
     *              left hand matrix.
     * 
     * @return the matrix product.
     */
    public abstract M multiply(M right);
    
    /**
     * Returns the element at {@code y}th row, {@code x}th column.
     * 
     * @param x the {@code X}-coordinate of the element.
     * @param y the {@code Y}-coordinate of the element.
     * 
     * @return the matrix element at specified coordinates.
     */
    public abstract E get(int x, int y);
    
    /**
     * Sets the value {@code value} at {@code y}th row, {@code x}th column.
     * 
     * @param x     the {@code X}-coordinate of the value.
     * @param y     the {@code Y}-coordinate of the value.
     * @param value the value to set.
     */
    public abstract void set(int x, int y, E value);
    
    /**
     * Checks whether {@code o} is an abstract matrix and has the same content
     * as this matrix.
     * 
     * @param o the matrix to check against.
     * @return {@code true} only if the two matrices are equal.
     */
    @Override
    public boolean equals(Object o) {
        AbstractMatrix<M, E> other = (AbstractMatrix<M, E>) o;
        
        if (width != other.width || height != other.height) {
            return false;
        }
        
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                if (!get(x, y).equals(other.get(x, y))) {
                    return false;
                }
            }
        }
        
        return true;
    }
    
    protected static int checkWidth(int widthCandidate) {
        if (widthCandidate == 0) {
            throw new IllegalArgumentException("Matrix width is zero.");
        }
        
        if (widthCandidate < 0) {
            throw new IllegalArgumentException(
                    "Matrix width is negative: " + widthCandidate);
        }
        
        return widthCandidate;
    }
    
    protected static int checkHeight(int heightCandidate) {
        if (heightCandidate == 0) {
            throw new IllegalArgumentException("Matrix width is zero.");
        }
        
        if (heightCandidate < 0) {
            throw new IllegalArgumentException(
                    "Matrix width is negative: " + heightCandidate);
        }
        
        return heightCandidate;
    }
    
    protected void checkCoordinates(int x, int y) {
        checkX(x);
        checkY(y);
    }
    
    private void checkX(int x) {
        if (x < 0) {
            throw new IndexOutOfBoundsException(
                    "X-coordinate is negative: " + x);
        }
        
        if (x >= width) {
            throw new IndexOutOfBoundsException(
                    "X-coordinate is too large: " 
                            + x
                            + ". Must be at most " 
                            + (width - 1) 
                            + ".");
        }
    }
    
    private void checkY(int y) {
        if (y < 0) {
            throw new IndexOutOfBoundsException(
                    "Y-coordinate is negative: " + y);
        }
        
        if (y >= height) {
            throw new IndexOutOfBoundsException(
                    "Y-coordinate is too large: " 
                            + y
                            + ". Must be at most " 
                            + (height - 1) 
                            + ".");
        }
    }
    
    protected void checkMatrixHaveSameDimensions(M matrix1, M matrix2) {
        if (matrix1.getWidth() != matrix2.getWidth()) {
            throw new MatricesNotAddableException(
                    "Matrix widths mismatch: " 
                            + matrix1.getWidth() 
                            + " vs " 
                            + matrix2.getWidth() 
                            + ".");
        }
        
        if (matrix1.getHeight() != matrix2.getHeight()) {
            throw new MatricesNotAddableException(
                    "Matrix heights mismatch: " 
                            + matrix1.getHeight() 
                            + " vs " 
                            + matrix2.getHeight() 
                            + ".");
        }
    }
        
    protected void checkMatricesCanBeMultiplied(M leftMatrix, M rightMatrix) {
        if (leftMatrix.getWidth() != rightMatrix.getHeight()) {
            throw new MatricesNotMultipliableException(
                    "Cannot multiply the matrices. Width of left matrix is "
                            + leftMatrix.getWidth() 
                            + ", the height of the right matrix is " 
                            + rightMatrix.getHeight() 
                            + ".");
        }
    }
}

com.github.coderodde.math.linear.matrix.DenseMatrix2D.java:

package com.github.coderodde.math.linear.matrix;

/**
 * This class implements a (dense) matrix stored as a two-dimensional array.
 * 
 * @param <E> the matrix element type.
 * @author Rodion "rodde" Efremov
 * @version 1.6 (Aug 13, 2023)
 * @since 1.6 (Aug 13, 2023)
 */
public class DenseMatrix2D<E> extends AbstractMatrix<DenseMatrix2D<E>, E> {
    
    /**
     * The actual matrix holding the elements.
     */
    private final E[][] data;

    /**
     * Constructs a new dense matrix that stores all the elements in a two-
     * dimensional array.
     * 
     * @param width         the width of this matrix.
     * @param height        the height of this matrix.
     * @param fieldElements the field element API object.
     */
    public DenseMatrix2D(int width, int height, FieldElement<E> fieldElements) {
        super(width, height, fieldElements);
        
        data = (E[][]) new Object[height][];
        
        for (int y = 0; y < height; y++) {
            data[y] = (E[]) new Object[width];
        }
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public E get(int x, int y) {
        checkCoordinates(x, y);
        return data[y][x];
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public void set(int x, int y, E value) {
        if (value == null || fieldElements.identity().equals(value)) {
            data[y][x] = fieldElements.identity(); 
        } else {
            data[y][x] = value;
        }
    }
    
    /**
     * {@inheritDoc }
     */
    @Override
    public void negate() {
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                set(x, y, fieldElements.negate(get(x, y)));
            }
        }
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public DenseMatrix2D<E> immutableNegate() {
        DenseMatrix2D<E> ret = new DenseMatrix2D<>(width,
                                                   height,
                                                   fieldElements);
        
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                ret.set(x, y, fieldElements.negate(get(x, y)));
            }
        }
        
        return ret;
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public void add(DenseMatrix2D<E> other) {
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                set(x, y, fieldElements.add(get(x, y), other.get(x, y)));
            }
        }
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public DenseMatrix2D<E> immutableAdd(DenseMatrix2D<E> other) {
        DenseMatrix2D<E> ret = new DenseMatrix2D<>(width, 
                                                   height,
                                                   fieldElements);
        
        for (int y = 0; y < height; y++) {
            for (int x = 0; x < width; x++) {
                ret.set(x, y, fieldElements.add(get(x, y), other.get(x, y)));
            }
        }
        
        return ret;
    }

    /**
     * {@inheritDoc }
     */
    @Override
    public DenseMatrix2D<E> multiply(DenseMatrix2D<E> right) {
        DenseMatrix2D<E> ret = new DenseMatrix2D<>(right.getWidth(),
                                                   this.getHeight(),
                                                   fieldElements);
        
        for (int row = 0; row < height; row++) {
            for (int col = 0; col < right.getWidth(); col++) {
                ret.set(col, row, combineRowColumn(col, row, right));
            }
        }
        
        return ret;
    }
    
    private E combineRowColumn(int col, int row, DenseMatrix2D<E> right) {
        E sum = fieldElements.identity();
        
        for (int x = 0; x < width; x++) {
            E product = fieldElements.multiply(get(x, row), right.get(col, x));
            sum = fieldElements.add(sum, product);
        }
        
        return sum;
    }
}

com.github.coderodde.math.linear.matrix.SparseMatrix.java:

package com.github.coderodde.math.linear.matrix;

import java.util.HashMap;
import java.util.Map;

/**
 * This class implements a sparse matrix.
 * 
 * @param <E> the matrix element type.
 * @author Rodion "rodde" Efremov
 * @version 1.6 (Aug 13, 2023)
 * @since 1.6 (Aug 13, 2023)
 */
public class SparseMatrix<E> extends AbstractMatrix<SparseMatrix<E>, E> {
    
    private final Map<Integer, Map<Integer, E>> dataXY = new HashMap<>();
    private final Map<Integer, Map<Integer, E>> dataYX = new HashMap<>();
    
    public SparseMatrix(int width, int height, FieldElement<E> fieldElements) {
        super(width, height, fieldElements);
    }

    @Override
    public void negate() {
        for (Map.Entry<Integer, Map<Integer, E>> entry1 : dataXY.entrySet()) {
            for (Map.Entry<Integer, E> entry2 : entry1.getValue().entrySet()) {
                entry2.setValue(fieldElements.negate(entry2.getValue()));
            }
        }
    }

    @Override
    public SparseMatrix<E> immutableNegate() {
        SparseMatrix<E> ret = new SparseMatrix<>(width, height, fieldElements);
        
        for (Map.Entry<Integer, Map<Integer, E>> entry1 : dataXY.entrySet()) {
            int x = entry1.getKey();
            
            for (Map.Entry<Integer, E> entry2 : entry1.getValue().entrySet()) {
                int y = entry2.getKey();
                
                ret.set(x, y, fieldElements.negate(entry2.getValue()));
            }
        }
        
        return ret;
    }

    @Override
    public void add(SparseMatrix<E> other) {
        for (Map.Entry<Integer, Map<Integer, E>> entry1
                : other.dataXY.entrySet()) {
            int x = entry1.getKey();
            
            for (Map.Entry<Integer, E> entry2 : entry1.getValue().entrySet()) {
                int y = entry2.getKey();
                
                set(x, y, fieldElements.add(get(x, y), entry2.getValue()));
            }
        }
    }

    @Override
    public SparseMatrix<E> immutableAdd(SparseMatrix<E> other) {
        SparseMatrix<E> ret = new SparseMatrix<>(width, height, fieldElements);
        
        for (Map.Entry<Integer, Map<Integer, E>> entry1 : dataXY.entrySet()) {
            int x = entry1.getKey();
            
            for (Map.Entry<Integer, E> entry2 : entry1.getValue().entrySet()) {
                int y = entry2.getKey();
                
                ret.set(x, y, entry2.getValue());
            }
        }
        
        for (Map.Entry<Integer, Map<Integer, E>> entry1
                : other.dataXY.entrySet()) {
            int x = entry1.getKey();
            
            for (Map.Entry<Integer, E> entry2 : entry1.getValue().entrySet()) {
                int y = entry2.getKey();
                
                ret.set(x, 
                        y, 
                        fieldElements.add(ret.get(x, y), 
                                          other.get(x, y)));
            }
        }
        
        return ret;
    }

    @Override
    public SparseMatrix<E> multiply(SparseMatrix<E> right) {
        checkMatricesCanBeMultiplied(this, right);
        SparseMatrix<E> ret = new SparseMatrix<>(width, height, fieldElements);
    
        for (int leftRow = 0; leftRow < height; leftRow++) {
            for (int rightColumn = 0; 
                    rightColumn < right.width;
                    rightColumn++) {
                
                ret.set(rightColumn,
                        leftRow, 
                        combineRowCol(dataYX.get(leftRow),
                                      right.dataXY.get(rightColumn)));
            }
        }
        
        return ret;
    }

    @Override
    public E get(int x, int y) {
        if (!dataXY.containsKey(x)) {
            return fieldElements.identity();
        }
        
        return dataXY.get(x).getOrDefault(y, fieldElements.identity());
    }

    @Override
    public void set(int x, int y, E value) {
        if (value == null || value.equals(fieldElements.identity())) {
            deleteZeroEntry(x, y);
        } else {
            updateEntry(x, y, value);
        }
    }
    
    private void deleteZeroEntry(int x, int y) {
        if (dataXY.containsKey(x)) {
            dataXY.get(x).remove(y);
            
            if (dataXY.get(x).isEmpty()) {
                dataXY.remove(x);
            }
        }
        
        if (dataYX.containsKey(y)) {
            dataYX.get(y).remove(x);
            
            if (dataYX.get(y).isEmpty()) {
                dataYX.remove(y);
            }
        }
    }
    
    private void updateEntry(int x, int y, E value) {
        if (!dataXY.containsKey(x)) {
            Map<Integer, E> subMap = new HashMap<>();
            subMap.put(y, value);
            dataXY.put(x, subMap);
        } else {
            dataXY.get(x).put(y, value);
        }
        
        if (!dataYX.containsKey(y)) {
            Map<Integer, E> subMap = new HashMap<>();
            subMap.put(x, value);
            dataYX.put(y, subMap);
        } else {
            dataYX.get(y).put(x, value);
        }
    }
    
    private E combineRowCol(Map<Integer, E> map1, Map<Integer, E> map2) {
        E sum = fieldElements.identity();
        
        if (map1 == null || map2 == null) {
            return sum;
        }
        
        if (map1.size() < map2.size()) {
            return combineRowCol(map2, map1);
        }
        
        for (Map.Entry<Integer, E> entry : map2.entrySet()) {
            int a = entry.getKey();
            E rowValue = entry.getValue();
            
            if (map1.containsKey(a)) {
                E columnValue = map1.get(a);
                E product = fieldElements.multiply(columnValue, rowValue);
                sum = fieldElements.add(sum, product);
            }
        }
        
        return sum;
    }
}

Typical demo output

The demonstration program outputs something like that:

Warming up...
Benchmarking...
Created dense matrix in 1 ms.
Created sparse matrix in 1 ms.
Dense matrix addition in 15 ms.
Sparse matrix addition in 2 ms.
Addition matches: true
Dense matrix multiplication in 15572 ms.
Sparse matrix multiplication in 185 ms.
Multiplication matches: true

So we see that sparse matrix multiplication is much faster than the dense matrix multiplication.

Critique request

Can I improve anything here? Please tell me anything that comes to mind.

Answer 1

Dense matrix multiplication in 15572 ms.

Going by the demo code I assume this was the product between two 1000x1000 matrices. By my estimate based on old code the raw operations by themselves could take somewhere around 70ms under good conditions (on a PC from a decade ago, single threaded), which sounds mad compared to 15572ms, but the conditions that enabled that are incompatible with this generic API.

• You probably cannot use SIMD. Java has SIMD APIs these days, but you're dealing with generic elements. Even if you cheat with instanceof and write some type-specific implementations, there would be no SIMD-friendly access to the underlying values. If there were type-specific matrix types (so that a matrix of integers has a raw int[] to work with) then it could be done, I'm not necessarily recommending that, you can think about it.
• You can use loop-tiling to improve the memory access pattern partially. The extra indirection throws a bit of a wrench into the memory access pattern anyway, but I expect some gain from that.

There are not many algorithms implemented (yet?) so not much to really do with these matrices. No transpose even. As for other algorithms I suppose a blocker is that:

public interface FieldElement<E> extends Negable<E>, 
                                         Addable<E>, 
                                         Multipliable<E>,
                                         Zero<E> {
 }

.. is not a field, just a ring. It can add, multiply, and negate, but not take multiplicative inverses. That is relatively annoying in linear algebra.

You could do more with actual fields, eg LU-decomposition and solving Ax=b. That includes not just float/double and complex, but also finite fields - for that aspect perhaps you can draw some inspiration from FFLAS-FFPACK

BTW here's some video (first in a series) about sparse matrix algorithms, maybe you can use it to draw inspiration from it too: https://www.youtube.com/watch?v=1dGRTOwBkQs

---

@Override
public SparseMatrix<E> multiply(SparseMatrix<E> right) {
    checkMatricesCanBeMultiplied(this, right);
    SparseMatrix<E> ret = new SparseMatrix<>(width, height, fieldElements);

    for (int leftRow = 0; leftRow < height; leftRow++) {
        for (int rightColumn = 0; 
                rightColumn < right.width;
                rightColumn++) {
            
            ret.set(rightColumn,
                    leftRow, 
                    combineRowCol(dataYX.get(leftRow),
                                  right.dataXY.get(rightColumn)));
        }
    }
    
    return ret;
}

Having dataXY and dataYX is a neat trick, enabling simple efficient sparse matrix multiplication. You don't actually need to store the matrix twice though, there are algorithms for multiplying sparse matrices that are represented in only column format or only row format too (or two matrices with opposite representation, but that's easy, then you get the same algorithm that you already have). There are various different approaches.

Here's one example, assuming both matrices are stored by column: take a column from the right matrix, for each non-zero element in it, take the corresponding column from the left matrix and generate all products. Add the products into their corresponding locations in the result matrix, creating new entries as necessary. This is normally an annoying approach (ie when directly manipulating a matrix in CSC format) but with a map of maps it's actually fine. I'm not sure if I got this algorithm 100% right, but the approach is something like this anyway.

Answer 2

just some minor thoughts...

naming

I really lave how you name the variables when you iterate over the matrix, seen in many ways in your code. I really apprechiate that it is NOT i / j but instead x / y to present the location in your matrix - great!

aside in one loop in DenseMatrix2D, where you switch from x / y to row / column - that breaks up the consistency a bit

using more ide support

you could benefit far more from the features of an IDE - I am not sure if you use one, but it would help in formatting, which is kind of "unexpected" in some case and not so consistent.

Aside of that you would have seen a double semi-colon in AbstractMatrix which is not a issue but a smell. this.fieldElements = Objects.requireNonNull(fieldElements);;

incomplete code

I cannot find code for MatricesNotAddableException or MatricesNotMultipliableException which is sad. and i am not sure if you really need a new Exception and do not rely on an well-know ArithmeticException - but that is very opinion based ^^

javadoc

/**
 * This class implements a sparse matrix.
 * 
 * @param <E> the matrix element type.
 * @author Rodion "rodde" Efremov
 * @version 1.6 (Aug 13, 2023)
 * @since 1.6 (Aug 13, 2023)
 */
public class SparseMatrix<E> extends AbstractMatrix<SparseMatrix<E>, E> {...}

could as well be skipped... on some others i am also not sure if they are not bloaters rather than helpers - again very opinion based