Sparse matrices / arrays in Java
I'm working on a project, written in Java, which requires that I build a very large 2-D sparse array. Very sparse, if that makes a difference. Anyway: the most crucial aspect for this application is efficency in terms of time (assume loads of memory, though not nearly so unlimited as to allow me to use a standard 2-D array -- the key range is in the billions in both dimensions).
Out of the kajillion cells in the array, there will be several hundred thousand cells which contain an object. I need to be able to modify cell contents VERY quickly.
Anyway: Does anyone know a particularly good library for this purpose? It would have to be Berkeley, LGPL or similar license (no GPL, as the product can't be entirely open-sourced). Or if there's just a very simple way to make a homebrew sparse array object, that'd be fine too.
I'm considering MTJ, but haven't heard any opinions on its quality.
Solution 1:
Sparsed arrays built with hashmaps are very inefficient for frequently read data. The most efficient implementations uses a Trie that allows access to a single vector where segments are distributed.
A Trie can compute if an element is present in the table by performing only read-only TWO array indexing to get the effective position where an element is stored, or to know if its absent from the underlying store.
It can also provide a default position in the backing store for the default value of the sparsed array, so that you don't need ANY test on the returned index, because the Trie guarantees that all possible source index will map at least to the default position in the backing store (where you'll frequently store a zero, or an empty string or a null object).
There exists implementations that support fast-updatable Tries, with an otional "compact()" operation to optimze the size of the backing store at end of multiple operations. Tries are MUCH faster than hashmaps, because they don't need any complex hashing function, and don't need to handle collisions for reads (With Hashmaps, you have collision BOTH for reading and for writing, this requires a loop to skip to the next candidate position, and a test on each of them to compare the effective source index...)
In addition, Java Hashmaps can only index on Objects, and creating an Integer object for each hashed source index (this object creation will be needed for every read, not just writes) is costly in terms of memory operations, as it stresses the garbage collector.
I really hoped that the JRE included an IntegerTrieMap<Object> as the default implementation for the slow HashMap<Integer, Object> or LongTrieMap<Object> as the default implementation for the even slower HashMap<Long, Object>... But this is still not the case.
You may wonder what is a Trie?
It's just a small array of integers (in a smaller range than the full range of coordinates for your matrix) that allows mapping the coordinates into an integer position in a vector.
For example suppose you want a 1024*1024 matrix containing only a few non zero values. Instead of storing that matrix in a array containing 1024*1024 elements (more than 1 million), you may just want to split it into subranges of size 16*16, and you'll just need 64*64 such subranges.
In that case, the Trie index will contain just 64*64 integers (4096), and there will be at least 16*16 data elements (containing the default zeroes, or the most common subrange found in your sparse matrix).
And the vector used to store the values will contain only 1 copy for subranges that are equal with each other (most of them being full of zeroes, they will be represented by the same subrange).
So instead of using a syntax like matrix[i][j], you'd use a syntax like:
trie.values[trie.subrangePositions[(i & ~15) + (j >> 4)] +
((i & 15) << 4) + (j & 15)]
which will be more conveniently handled using an access method for the trie object.
Here is an example, built into a commented class (I hope it compiles OK, as it was simplified; signal me if there are errors to correct):
/**
* Implement a sparse matrix. Currently limited to a static size
* (<code>SIZE_I</code>, <code>SIZE_I</code>).
*/
public class DoubleTrie {
/* Matrix logical options */
public static final int SIZE_I = 1024;
public static final int SIZE_J = 1024;
public static final double DEFAULT_VALUE = 0.0;
/* Internal splitting options */
private static final int SUBRANGEBITS_I = 4;
private static final int SUBRANGEBITS_J = 4;
/* Internal derived splitting constants */
private static final int SUBRANGE_I =
1 << SUBRANGEBITS_I;
private static final int SUBRANGE_J =
1 << SUBRANGEBITS_J;
private static final int SUBRANGEMASK_I =
SUBRANGE_I - 1;
private static final int SUBRANGEMASK_J =
SUBRANGE_J - 1;
private static final int SUBRANGE_POSITIONS =
SUBRANGE_I * SUBRANGE_J;
/* Internal derived default values for constructors */
private static final int SUBRANGES_I =
(SIZE_I + SUBRANGE_I - 1) / SUBRANGE_I;
private static final int SUBRANGES_J =
(SIZE_J + SUBRANGE_J - 1) / SUBRANGE_J;
private static final int SUBRANGES =
SUBRANGES_I * SUBRANGES_J;
private static final int DEFAULT_POSITIONS[] =
new int[SUBRANGES](0);
private static final double DEFAULT_VALUES[] =
new double[SUBRANGE_POSITIONS](DEFAULT_VALUE);
/* Internal fast computations of the splitting subrange and offset. */
private static final int subrangeOf(
final int i, final int j) {
return (i >> SUBRANGEBITS_I) * SUBRANGE_J +
(j >> SUBRANGEBITS_J);
}
private static final int positionOffsetOf(
final int i, final int j) {
return (i & SUBRANGEMASK_I) * MAX_J +
(j & SUBRANGEMASK_J);
}
/**
* Utility missing in java.lang.System for arrays of comparable
* component types, including all native types like double here.
*/
public static final int arraycompare(
final double[] values1, final int position1,
final double[] values2, final int position2,
final int length) {
if (position1 >= 0 && position2 >= 0 && length >= 0) {
while (length-- > 0) {
double value1, value2;
if ((value1 = values1[position1 + length]) !=
(value2 = values2[position2 + length])) {
/* Note: NaN values are different from everything including
* all Nan values; they are are also neigher lower than nor
* greater than everything including NaN. Note that the two
* infinite values, as well as denormal values, are exactly
* ordered and comparable with <, <=, ==, >=, >=, !=. Note
* that in comments below, infinite is considered "defined".
*/
if (value1 < value2)
return -1; /* defined < defined. */
if (value1 > value2)
return 1; /* defined > defined. */
if (value1 == value2)
return 0; /* defined == defined. */
/* One or both are NaN. */
if (value1 == value1) /* Is not a NaN? */
return -1; /* defined < NaN. */
if (value2 == value2) /* Is not a NaN? */
return 1; /* NaN > defined. */
/* Otherwise, both are NaN: check their precise bits in
* range 0x7FF0000000000001L..0x7FFFFFFFFFFFFFFFL
* including the canonical 0x7FF8000000000000L, or in
* range 0xFFF0000000000001L..0xFFFFFFFFFFFFFFFFL.
* Needed for sort stability only (NaNs are otherwise
* unordered).
*/
long raw1, raw2;
if ((raw1 = Double.doubleToRawLongBits(value1)) !=
(raw2 = Double.doubleToRawLongBits(value2)))
return raw1 < raw2 ? -1 : 1;
/* Otherwise the NaN are strictly equal, continue. */
}
}
return 0;
}
throw new ArrayIndexOutOfBoundsException(
"The positions and length can't be negative");
}
/**
* Utility shortcut for comparing ranges in the same array.
*/
public static final int arraycompare(
final double[] values,
final int position1, final int position2,
final int length) {
return arraycompare(values, position1, values, position2, length);
}
/**
* Utility missing in java.lang.System for arrays of equalizable
* component types, including all native types like double here.
*/
public static final boolean arrayequals(
final double[] values1, final int position1,
final double[] values2, final int position2,
final int length) {
return arraycompare(values1, position1, values2, position2, length) ==
0;
}
/**
* Utility shortcut for identifying ranges in the same array.
*/
public static final boolean arrayequals(
final double[] values,
final int position1, final int position2,
final int length) {
return arrayequals(values, position1, values, position2, length);
}
/**
* Utility shortcut for copying ranges in the same array.
*/
public static final void arraycopy(
final double[] values,
final int srcPosition, final int dstPosition,
final int length) {
arraycopy(values, srcPosition, values, dstPosition, length);
}
/**
* Utility shortcut for resizing an array, preserving values at start.
*/
public static final double[] arraysetlength(
double[] values,
final int newLength) {
final int oldLength =
values.length < newLength ? values.length : newLength;
System.arraycopy(values, 0, values = new double[newLength], 0,
oldLength);
return values;
}
/* Internal instance members. */
private double values[];
private int subrangePositions[];
private bool isSharedValues;
private bool isSharedSubrangePositions;
/* Internal method. */
private final reset(
final double[] values,
final int[] subrangePositions) {
this.isSharedValues =
(this.values = values) == DEFAULT_VALUES;
this.isSharedsubrangePositions =
(this.subrangePositions = subrangePositions) ==
DEFAULT_POSITIONS;
}
/**
* Reset the matrix to fill it with the same initial value.
*
* @param initialValue The value to set in all cell positions.
*/
public reset(final double initialValue = DEFAULT_VALUE) {
reset(
(initialValue == DEFAULT_VALUE) ? DEFAULT_VALUES :
new double[SUBRANGE_POSITIONS](initialValue),
DEFAULT_POSITIONS);
}
/**
* Default constructor, using single default value.
*
* @param initialValue Alternate default value to initialize all
* positions in the matrix.
*/
public DoubleTrie(final double initialValue = DEFAULT_VALUE) {
this.reset(initialValue);
}
/**
* This is a useful preinitialized instance containing the
* DEFAULT_VALUE in all cells.
*/
public static DoubleTrie DEFAULT_INSTANCE = new DoubleTrie();
/**
* Copy constructor. Note that the source trie may be immutable
* or not; but this constructor will create a new mutable trie
* even if the new trie initially shares some storage with its
* source when that source also uses shared storage.
*/
public DoubleTrie(final DoubleTrie source) {
this.values = (this.isSharedValues =
source.isSharedValues) ?
source.values :
source.values.clone();
this.subrangePositions = (this.isSharedSubrangePositions =
source.isSharedSubrangePositions) ?
source.subrangePositions :
source.subrangePositions.clone());
}
/**
* Fast indexed getter.
*
* @param i Row of position to set in the matrix.
* @param j Column of position to set in the matrix.
* @return The value stored in matrix at that position.
*/
public double getAt(final int i, final int j) {
return values[subrangePositions[subrangeOf(i, j)] +
positionOffsetOf(i, j)];
}
/**
* Fast indexed setter.
*
* @param i Row of position to set in the sparsed matrix.
* @param j Column of position to set in the sparsed matrix.
* @param value The value to set at this position.
* @return The passed value.
* Note: this does not compact the sparsed matric after setting.
* @see compact(void)
*/
public double setAt(final int i, final int i, final double value) {
final int subrange = subrangeOf(i, j);
final int positionOffset = positionOffsetOf(i, j);
// Fast check to see if the assignment will change something.
int subrangePosition, valuePosition;
if (Double.compare(
values[valuePosition =
(subrangePosition = subrangePositions[subrange]) +
positionOffset],
value) != 0) {
/* So we'll need to perform an effective assignment in values.
* Check if the current subrange to assign is shared of not.
* Note that we also include the DEFAULT_VALUES which may be
* shared by several other (not tested) trie instances,
* including those instanciated by the copy contructor. */
if (isSharedValues) {
values = values.clone();
isSharedValues = false;
}
/* Scan all other subranges to check if the position in values
* to assign is shared by another subrange. */
for (int otherSubrange = subrangePositions.length;
--otherSubrange >= 0; ) {
if (otherSubrange != subrange)
continue; /* Ignore the target subrange. */
/* Note: the following test of range is safe with future
* interleaving of common subranges (TODO in compact()),
* even though, for now, subranges are sharing positions
* only between their common start and end position, so we
* could as well only perform the simpler test <code>
* (otherSubrangePosition == subrangePosition)</code>,
* instead of testing the two bounds of the positions
* interval of the other subrange. */
int otherSubrangePosition;
if ((otherSubrangePosition =
subrangePositions[otherSubrange]) >=
valuePosition &&
otherSubrangePosition + SUBRANGE_POSITIONS <
valuePosition) {
/* The target position is shared by some other
* subrange, we need to make it unique by cloning the
* subrange to a larger values vector, copying all the
* current subrange values at end of the new vector,
* before assigning the new value. This will require
* changing the position of the current subrange, but
* before doing that, we first need to check if the
* subrangePositions array itself is also shared
* between instances (including the DEFAULT_POSITIONS
* that should be preserved, and possible arrays
* shared by an external factory contructor whose
* source trie was declared immutable in a derived
* class). */
if (isSharedSubrangePositions) {
subrangePositions = subrangePositions.clone();
isSharedSubrangePositions = false;
}
/* TODO: no attempt is made to allocate less than a
* fully independant subrange, using possible
* interleaving: this would require scanning all
* other existing values to find a match for the
* modified subrange of values; but this could
* potentially leave positions (in the current subrange
* of values) unreferenced by any subrange, after the
* change of position for the current subrange. This
* scanning could be prohibitively long for each
* assignement, and for now it's assumed that compact()
* will be used later, after those assignements. */
values = setlengh(
values,
(subrangePositions[subrange] =
subrangePositions = values.length) +
SUBRANGE_POSITIONS);
valuePosition = subrangePositions + positionOffset;
break;
}
}
/* Now perform the effective assignment of the value. */
values[valuePosition] = value;
}
}
return value;
}
/**
* Compact the storage of common subranges.
* TODO: This is a simple implementation without interleaving, which
* would offer a better data compression. However, interleaving with its
* O(N²) complexity where N is the total length of values, should
* be attempted only after this basic compression whose complexity is
* O(n²) with n being SUBRANGE_POSITIIONS times smaller than N.
*/
public void compact() {
final int oldValuesLength = values.length;
int newValuesLength = 0;
for (int oldPosition = 0;
oldPosition < oldValuesLength;
oldPosition += SUBRANGE_POSITIONS) {
int oldPosition = positions[subrange];
bool commonSubrange = false;
/* Scan values for possible common subranges. */
for (int newPosition = newValuesLength;
(newPosition -= SUBRANGE_POSITIONS) >= 0; )
if (arrayequals(values, newPosition, oldPosition,
SUBRANGE_POSITIONS)) {
commonSubrange = true;
/* Update the subrangePositions|] with all matching
* positions from oldPosition to newPosition. There may
* be several index to change, if the trie has already
* been compacted() before, and later reassigned. */
for (subrange = subrangePositions.length;
--subrange >= 0; )
if (subrangePositions[subrange] == oldPosition)
subrangePositions[subrange] = newPosition;
break;
}
if (!commonSubrange) {
/* Move down the non-common values, if some previous
* subranges have been compressed when they were common.
*/
if (!commonSubrange && oldPosition != newValuesLength) {
arraycopy(values, oldPosition, newValuesLength,
SUBRANGE_POSITIONS);
/* Advance compressed values to preserve these new ones. */
newValuesLength += SUBRANGE_POSITIONS;
}
}
}
/* Check the number of compressed values. */
if (newValuesLength < oldValuesLength) {
values = values.arraysetlength(newValuesLength);
isSharedValues = false;
}
}
}
Note: this code is not complete because it handles a single matrix size, and its compactor is limited to detect only common subranges, without interleaving them.
Also, the code does not determine where it is the best width or height to use for splitting the matrix into subranges (for x or y coordinates), according to the matrix size. It just uses the same static subrange sizes of 16 (for both coordinates), but it could be conveniently any other small power of 2 (but a non power of 2 would slow down the int indexOf(int, int) and int offsetOf(int, int) internal methods), independantly for both coordinates, and up to the maximum width or height of the matrix.ideally the compact() method should be able to determine the best fitting sizes.
If these splitting subranges sizes can vary, then there will be a need to add instance members for these subrange sizes instead of the static SUBRANGE_POSITIONS, and to make the static methods int subrangeOf(int i, int j) and int positionOffsetOf(int i, int j) into non static; and the initialization arrays DEFAULT_POSITIONSand DEFAULT_VALUES will need to be dropped or redefined differently.
If you want to support interleaving, basically you'll start by dividing the existing values in two of about the same size (both being a multiple of the minimum subrange size, the first subset possibly having one more subrange than the second one), and you'll scan the larger one at all successive positions to find a matching interleaving; then you'll try to match these values. Then you'll loop recursely by dividing the subsets in halves (also a multiple of the minimum subrange size) and you'll scan again to match these subsets (this will multiply the number of subsets by 2: you have to wonder if the doubled size of the subrangePositions index is worth the value compared to the existing size of the values to see if it offers an effective compression (if not, you stop there: you have found the optimum subrange size directly from the interleaving compression process). In that case; the subrange size will be mutable, during compaction.
But this code shows how you assign non-zero values and reallocate the data array for additional (non-zero) subranges, and then how you can optimize (with compact() after assignments have been performed using the setAt(int i, int j, double value) method) the storage of this data when there are duplicate subranges that may be unified within the data, and reindexed at the same position in the subrangePositions array.
Anyway, all the principles of a trie are implemented there:
It is always faster (and more compact in memory, meaning better locality) to represent a matrix using a single vector instead of a double-indexed array of arrays (each one allocated separately). The improvement is visible in the
double getAt(int, int)method !You save a lot of space, but when assigning values it may take time to reallocate new subranges. For this reason, the subranges should not be too small or reallocations will occur too frequently for setting up your matrix.
It is possible to transform an initial large matrix automatically into a more compact matrix by detecting common subranges. A typical implementation will then contain a method such as
compact()above. However, if get() access is very fast and set() is quite fast, compact() may be very slow if there are lots of common subranges to compress (for example when substracting a large non-sparse randomly-filled matrix with itself, or multiplying it by zero: it will be simpler and much faster in that case to reset the trie by instanciating a new one and dropping the old one).Common subranges use common storage in the data, so this shared data must be read-only. If you must change a single value without changing the rest of the matrix, you must first make sure that it is referenced only one time in the
subrangePositionsindex. Otherwise you'll need to allocate a new subrange anywhere (conveniently at end) of thevaluesvector, and then store the position of this new subrange into thesubrangePositionsindex.
Note that the generic Colt library, though very good as it is, is not as good when working on sparse matrice, because it uses hashing (or row-compresssed) technics which do not implement the support for tries for now, despite it is an excellent optimization, which is both space-saving and time-saving, notably for the most frequent getAt() operations.
Even the setAt() operation described here for tries saves lot of time (the way is is implemented here, i.e. without automatic compaction after setting, which could still be implemented based on demand and estimated time where compaction would still save lot of storage space at the price of time): the time saving is proportional to the number of cells in subranges, and space saving is inversely proportional to the number of cells per subrange. A good compromize if then to use a subrange size such the number of cells per subrange is the square root of the total number of cells in a 2D matrix (it would be a cubic root when working with 3D matrice).
Hashing technics used in Colt sparse matrix implementations have the inconvenience that they add a lot of storage overhead, and slow access time due to possible collisions. Tries can avoid all collisions, and can then warranty to save linear O(n) time to O(1) time in the worst cases, where (n) is the number of possible collisions (which, in case of sparse matrix, may be up to the number of non-default-value cells in the matrix, i.e. up to the total number of size of the matrix times a factor proportional to the hashing filling factor, for a non-sparse i.e. full matrix).
The RC (row-compressed) technics used in Colt are nearer from Tries, but this is at another price, here the compression technics used, that have very slow access time for the most frequent read-only get() operations, and very slow compression for setAt() operations. In addition, the compression used is not orthogonal, unlike in this presentation of Tries where orthogonality is preserved. Tries would also be preserve this orthogonality for related viewing operations such as striding, transposition (viewed as a striding operation based on integer cyclic modular operations), subranging (and subselections in general, including with sorting views).
I just hope that Colt will be updated in some future to implement another implementation using Tries (i.e. TrieSparseMatrix instead of just HashSparseMatrix and RCSparseMatrix). The ideas are in this article.
The Trove implementation (based in int->int maps) are also based on hashing technics similar to the Colt's HashedSparseMatrix, i.e. they have the same inconvenience. Tries will be a lot faster, with a moderate additional space consumed (but this space can be optimized and become even better than Trove and Colt, in a deferred time, using a final compact()ion operation on the resulting matrix/trie).
Note: this Trie implementation is bound to a specific native type (here double). This is voluntary, because generic implementation using boxing types have a huge space overhead (and are much slower in acccess time). Here it just uses native unidimensional arrays of double rather than generic Vector. But it is certainly possible to derive a generic implementation as well for Tries... Unfortunately, Java still does not allow writing really generic classes with all the benefits of native types, except by writing multiple implementations (for a generic Object type or for each native type), and serving all these operation via a type factory. The language should be able to automatically instanciate the native implementations and build the factory automatically (for now it is not the case even in Java 7, and this is something where .Net still maintains its advantage for really generic types that are as fast as native types).
Solution 2:
Following framework to test Java Matrix Libraries, provides also a good list of these! https://lessthanoptimal.github.io/Java-Matrix-Benchmark/
Tested Libraries:
* Colt
* Commons Math
* Efficient Java Matrix Library (EJML)
* Jama
* jblas
* JScience (Older benchmarks only)
* Matrix Toolkit Java (MTJ)
* OjAlgo
* Parallel Colt
* Universal Java Matrix Package (UJMP)
Solution 3:
This seems to be simple.
You could use a binary tree of the data using row*maxcolums+column as an index.
To find item, you simply calculate row*maxcolums+column and binary search the tree looking for it, if it's not there, you can return null (it's О(log n) where n is the number of cells which contain an object).
Solution 4:
Not the fastest runtime solution probably, but the fastest I could come up with that seems to work. Create an Index class and use it as a key for a SortedMap, like:
SortedMap<Index, Object> entries = new TreeMap<Index, Object>();
entries.put(new Index(1, 4), "1-4");
entries.put(new Index(5555555555l, 767777777777l), "5555555555l-767777777777l");
System.out.println(entries.size());
System.out.println(entries.get(new Index(1, 4)));
System.out.println(entries.get(new Index(5555555555l, 767777777777l)));
My Index class looks like this (with some help from Eclipse code generator).
public static class Index implements Comparable<Index>
{
private long x;
private long y;
public Index(long x, long y)
{
super();
this.x = x;
this.y = y;
}
public int compareTo(Index index)
{
long ix = index.x;
if (ix == x)
{
long iy = index.y;
if (iy == y)
{
return 0;
}
else if (iy < y)
{
return -1;
}
else
{
return 1;
}
}
else if (ix < x)
{
return -1;
}
else
{
return 1;
}
}
public int hashCode()
{
final int PRIME = 31;
int result = 1;
result = PRIME * result + (int) (x ^ (x >>> 32));
result = PRIME * result + (int) (y ^ (y >>> 32));
return result;
}
public boolean equals(Object obj)
{
if (this == obj)
return true;
if (obj == null)
return false;
if (getClass() != obj.getClass())
return false;
final Index other = (Index) obj;
if (x != other.x)
return false;
if (y != other.y)
return false;
return true;
}
public long getX()
{
return x;
}
public long getY()
{
return y;
}
}
Solution 5:
You migth look at la4j (Linear Algebra for Java) library. It supports CRS (Compressed Row Storage) as well as CCS (Compressed Column Storage) internal representaions for sparse matrices. So, these are the most efficient and fast internal stuctures for sparse data.
Here is the brief example of using sparse matrices in la4j:
Matrix a = new CRSMatrix(new double[][]{ // 'a' - CRS sparse matrix
{ 1.0, 0.0, 3.0 },
{ 0.0, 5.0, 0.0 },
{ 7.0, 0.0. 9.0 }
});
Matrix b = a.transpose(); // 'b' - CRS sparse matrix
Matrix c = b.multiply(a, Matrices.CCS_FACTORY); // 'c' = 'b' * 'a';
// 'c' - CCS sparse matrix