## You are here

Submitted by h2b on 30. December 2015 - 21:00

# LinAlg -- The Linear-Algebra Scala Library

LinAlg provides data types and operations for algebraic vectors and matrices.

Vector and Matrix elements can be of arbitrary type, provided that a factory for that element type is available. Currently, factories exist for `Double`, `Float`, `Long`, `Int`, `Short`, `Byte` and `Char`.

Vector and matrix (row and column) indices can be any integer (to state more precisely, an index must be in the interval `[Mindex, Maxdex]`, as defined in the `Index` object). But only elements corresponding to a subset of that whole integer range actually are stored (the concrete elements), while all other elements by definition are zero (virtual elements).

This concept (which I got to value a long time ago with the ALGOL68 library prelude TORRIX [1]) does not only provide a natural approach to operations with vectors and matrixes of different index ranges, but also prevents from index-out-of-bounds exceptions.

Concrete index ranges (those indices that correspond to concrete elements) implicitly or explicitly are assigned when a vector or matrix is created. By default, it starts at 1 and extends to the number of elements specified, but this is customizable.

`Vector`s and `Matrix`es both are immutable by design, so there exist no `update` methods or the like (you cannot do `v(i) = something`). There are, however, builder classes that allow you to build `Vector`s and `Matrix`es element by element.

Currently, operations on vectors and matrices require identical element types. For instance, you can add a `Vector[Int]` to another `Vector[Int]`, but you cannot add a `Vector[Int]` to another `Vector[Double]`.

## L1 Indices

### L1-1 Vector

All vectors have an index range from `Index.Mindex` to `Index.Maxdex` (both inclusive). In addition, each individual vector `v` has a concrete index range from `v.index.low` to `v.index.high` (both inclusive).

``` Mindex       low    high        Maxdex

0   0   0   x   x   x   0   0   0

virtual     concrete      virtual
```

Only the values corresponding to the concrete indices are stored while the rest are virtual zeroes. The representation of a virtual zero depends on the element type (e.g., it is `0.0` for `Double` or `0` for `Int`). Operations on vectors treat virtual elements like concrete ones, i.e., as if they were stored as concrete values.

The values of `Mindex` and `Maxdex` are implementation-dependent, but they will be very close to `Int.MinValue` and `Int.MaxValue`, respectively. So, as long as you don't explore the very edges of the whole `Int` range, you don't have to bother about index-out-of-bounds exceptions.

### L1-2 Matrix

All matrices have both a row and column index range from `Index.Mindex` to `Index.Maxdex` (both inclusive). In addition, each individual matrix `a` has a concrete row index range from `a.index.dim1.low` to `a.index.dim1.high` and a concrete column index range from `a.index.dim2.low` to `a.index.dim2.high` (all inclusive).

``` Mindex    dim2.low dim2.high    Maxdex

0   0   0   0   0   0   0   0   0   Mindex

0  vir  0   x   x   x   0  vir  0   dim1.low

0  tu   0   concrete    0   tu  0

0  al   0   x   x   x   0   al  0   dim1.high

0   0   0   0   0   0   0   0   0   Maxdex

Mindex    dim2.low dim2.high    Maxdex
```

Only the values corresponding to the concrete indices are stored while the rest are virtual zeroes. The representation of a virtual zero depends on the element type (e.g., it is `0.0` for `Double` or `0` for `Int`). Depending on the implementation it is even likely that more zeroes as shown actually are virtual and not stored. Operations on matrices treat virtual elements like concrete ones, i.e., as if they were stored as concrete values.

The values of `Mindex` and `Maxdex` are implementation-dependent, but they will be very close to `Int.MinValue` and `Int.MaxValue`, respectively. So, as long as you don't explore the very edges of the whole `Int` range, you don't have to bother about index-out-of-bounds exceptions.

## L2 Creation

### L2-1 Vector

The most common way to create a vector is by simply enumerating its elements. You get a vector according to the type of the element sequence with a concrete index range starting at 1:

``````val vi = Vector(1,2,3) //vi is Vector[Int] with v(1)==1, v(2)==2, v(3)==3 and v(i)==0 for all other i in [Mindex, Maxdex]
val vd = Vector(1.0,2.0,3.0) //vd is Vector[Double] with v(1)==1.0, v(2)==2.0, v(3)==3.0 and v(i)==0.0 for all other i in [Mindex, Maxdex]``````

To get a vector with customized index range, use `Vector.at` and provide `index.low` as a separate (first) parameter:

``val vat0 = Vector.at(0)(1,2,3) //v0 is Vector[Int] with v(0)==1, v(1)==2, v(2)==3 and v(i)==0 for all other i in [Mindex, Maxdex]``

Or provide a generating function `f: Int => E` (that maps indices to elements of type `E`) with explicit concrete index bounds:

``val vf = Vector((i: Int) => i*i: Double, 1, 3) //vf is Vector[Double] with v(1)==1.0, v(2)==4.0, v(3)==9.0 and v(i)==0.0 for all other i in [Mindex, Maxdex]``

If you cannot enumerate the vector elements in advance use a vector builder:

``````val vbuilder1 = Vector.newBuilder[Double]
vbuilder1 += 1.0
vbuilder1 += 2.0
vbuilder1 += 3.0
val vb1 = vbuilder1.result()``````

or

``````val vbuilder2 = Vector.newBuilder[Double]
vbuilder2(1) = 1.0
vbuilder2(3) = 3.0
val vb2 = vbuilder2.result()``````

In the first case you get a `Vector[Double](1.0,2.0,3.0)`, in the second case the result is `Vector[Double](1.0,0.0,3.0)`. Note that the second vector implicitly got a (concrete) element `vb2(2)==0.0`.

### L2-2 Matrix

Matrices are built from rows which are established by vectors.

Create a matrix by enumerating its row vectors. The matrix' type is derived from the given vectors, its concrete row index starts at 1 by default and the concrete column range is implied from the row vectors again:

``val a = Matrix(Vector(11,12,13), Vector(21,22,23)) //a is Matrix[Int] with a(1,1)==11, a(1,2)==12, ..., a(2,3)==23 and a(i,j)==0 for all not explicit or implicit mentioned (i,j) in [Mindex, Maxdex]``

Customized vector indices will be honored as in

``val b = Matrix(Vector.at(0)(10,11,12), Vector.at(0)(20,21,22)) //b is Matrix[Int] with b(1,0)==10, b(1,1)==11, ..., b(2,2)==22 and b(i,j)==0 for all not explicit or implicit mentioned (i,j) in [Mindex, Maxdex]``

To customize the row range use `Matrix.atRow` with specified low index:

``val c = Matrix.atRow(0)(Vector(1,2,3), Vector(11,12,13)) //b is Matrix[Int] with a(0,1)==1, a(0,2)==2, ..., a(1,3)==13 and a(i,j)==0 for all not explicit or implicit mentioned (i,j) in [Mindex, Maxdex]``

Or provide a generating function `f: Int => Vector[E]` (that maps indices to elements of type `Vector[E]`, i.e., row vectors) with explicit concrete row-index bounds:

``val af = Matrix((i: Int) => Vector(i*10+1,i*10+2,i*10+3), 1, 2) //af is equal to a, see above``

If you cannot enumerate the matrix rows in advance use a matrix builder:

``````val row1 = Vector(1,2,3)
val row2 = Vector(4,5,6)
val mbuilder1 = Matrix.newBuilder[Int]
mbuilder1 += row1
mbuilder1 += row2
val ab1 = mbuilder1.result()``````

or

``````val mbuilder2 = Matrix.newBuilder[Int]
mbuilder2(1) = row1
mbuilder2(3) = row2
val ab2 = mbuilder2.result()``````

In the first case you get a `Matrix[Int](Vector(1,2,3),Vector(4,5,6))`, in the second case the result is `Matrix[Int](Vector(1,2,3),Vector(),Vector(4,5,6)`. Note that the second matrix implicitly got a (concrete but empty) vector for the second row.

## L3 Operations

All operations on vectors and matrices require an identical element type.

Beyond, due to the concept of virtual elements (see above), operations on vectors or matrices with different index ranges resolve properly.

### L3-1 Vector

Let `u` and `v` be vectors of some type `Vector[E]`, let `s` be a scalar of type `E` and `i` be an integer, then the following operations can be done:

``````u.length //size of the concrete index range
u.isZero //if u contains no elements other than zero
u(i)     //element with index i, zero if i outside concrete index range
+u       //unary plus
-u       //unary minus
u + v    //sum
u - v    //difference
u * s    //multiplication by scalar
s * u    //multiplication by scalar (*)
u * v    //scalar product (result of type `E`)
u.norm   //Euklidian norm (result always of type `Double`)
u @@ i   //copy of u with lower index bound i``````

(*) To use this, you must import the package object by `import de.h2b.scala.lib.math.linalg._`. Beware that -- in contrast to the other operators -- this implementation of `s * u` always returns a `Vector[E]`; thus, subtypes should bring along their own implementation of this operator.

### L3-2 Matrix

Let `a` and `b` be matrices of some type `Matrix[E]`, let `v` be a vector of type `Vector[E]`, let `s` be a scalar of type `E` and `i` and `j` be integer, then the following operations can be done:

``````a.height    //number of concrete rows
a.width     //number of concrete columns
a.isSquare  //if a is a square matrix (regarding to the concrete index range)
a.isZero    //if a contains no elements other than zero
a(i,j)      //element with index (i,j), zero if (i,j) outside concrete index range
a(i)        //same as a.row(i)
a.row(i)    //the row vector with index i, zero vector if i outside concrete index range
a.col(j)    //the column vector with index j, zero vector if j outside concrete index range
+a          //unary plus
-a          //unary minus
a + b       //sum
a - b       //difference
a * b       //product
a * v       //product
v ** a      //product (*) (**)
a * s       //multiplication by scalar
s * a       //multiplication by scalar (*)
a.transpose //transpose matrix
a.rowSum    //vector of row sums of a
a.colSum    //vector of column sums of a
a.atRow(i)  //copy of a with lower row-index bound i
a.atCol(j)  //copy of a with lower column-index bound j (all row vectors are shifted to this bound)
a @@ (i,j)  //copy of a with lower row-index bound i and lower column-index bound j (all row vectors are shifted to j)``````

(*) To use this, you must import the package object by `import de.h2b.scala.lib.math.linalg._`. Beware that -- in contrast to the other operators -- these implementations of `s * a` and `v ** a` always return a `Matrix[E]` and `Vector[E]`, respectively; thus, subtypes should bring along their own implementations of these operators.

(**) This operator should be named *, but there is an ambiguity with scalar ops. Until this problem is solved, ** must be used.

## L4 Iterators

Both vector and matrix extend `Iterable` and thus provide iterators over its elements. In addition, you can use all methods inherited from `Iterable` like `map` or `foreach`.

### L4-1 Vector

`Vector[E]` extends `Iterable[E]` and provides an iterator over its elements:

``for (x <- Vector(1,2,3)) print(x) //prints 123``

### L4-2 Matrix

`Matrix[E]` extends `Iterable[Vector[E]]`.

A matrix `a` provides three iterators over its elements:

``````a.rowIterator  //iterates over the row vectors of a
a.colIterator  //iterates over the column vectors of a
a.elemIterator //iterates over the scalar elements of a``````

The first two yield vectors over which you can iterate (`a(i,)` or `a(,j)`) while the third one iterates over all bare elements (`a(i,j)`).

The `rowIterator` is the default one, thus

``````for (x <- Matrix(Vector(11,12,13), Vector(21,22,23))) println(x)
//prints
//(11,12,13)@1
//(21,22,23)@1``````

## L5 Sparse Vectors and Matrices

Since version 3.0.0 sparse vectors and matrices are provided. They come with storage engines that only store values different from zero. In addition, they are integrated seamless into the LinAlg system, i.e., they are `Vector`s and `Matrix`es and inherit all operations from them. Also, operations between sparse and ordinary vectors and matrices can be intermixed.

### L5-1 Vector

A sparse vector is created by a sequence of index-value mappings:

``val sv = Vector(-3 → -3, 1 → 1, 3 → 3)``

### L5-2 Matrix

A sparse matrix is created by a sequence of index-value mappings for its row vectors:

``val sm = Matrix(-1 → Vector(-11,-12,-13), +1 → Vector(11,12,13))``

Note that such a matrix is sparse regarding to its rows but not necessarily to its columns. Of course, the row vectors itself can be (but not need to be) sparse vectors.

## L6 Other Features

### L6-1 Equality

Two vectors or matrices `x` and `y` are considered to be equal (i.e., `x==y`) if and only if

1. `x` and `y` are of the same type (which includes the element type),
2. `x` and `y` have the same concrete index ranges and
3. `x` and `y` have the same elements at equal index positions.

This means, e.g., `Vector(1,2) != Vector(1,2,0)` (because of rule 2) and of course `Vector(1,2) != Vector(1,3)` (because of rule 3). With these examples, only something like `Vector(1,2) == Vector(1,2)` is true. The same applies to matrices.

For both vectors and matrices, there exists a similarity operator `~~` that lowers these restrictions and disregards rule 2 above. Two vectors or matrices `x` and `y` are considered to be similar (i.e., `x~~y`) if and only if

1. `x` and `y` are of the same type (which includes the element type) and
2. `x` and `y` have the same elements at equal index positions.

This means, e.g., `Vector(1,2) ~~ Vector(1,2,0)`. Of course, `! (Vector(1,2) ~~ Vector(1,3))` still holds and also `Vector(1,2) ~~ Vector(1,2)` is true. Again, the same applies to matrices.

### L6-2 Zero Vector and Matrix

To create zero vectors or matrices (i.e., those containing no concrete elements) just provide empty element lists but specify the element type:

``````val v0 = Vector[Double]()
val a0 = Matrix[Double]()``````

### L6-3 Subtyping

When subtyping the vector or matrix trait, of course all operations are inherited. Since version 2.0.0, a few measures ensure that these operations (as well as the ones inherited from the iterable trait which is extended by both the vector and matrix trait) will return proper result types.

For instance, let's implement a `Point` type by extending the vector trait. Our `Point` class is constructed from two double elements `x` and `y` and extends the `DoubleVector` class which is an implementation of `Vector[Double]`. For the result type of operations we mix in the `IterableLike` and `VectorLike` traits both with type parameters `Double` (the element type) and `Point` (the result type):

``````class Point private (val x: Double, val y: Double) extends
DoubleVector(1, Seq(x,y)) with IterableLike[Double, Point] with
VectorLike[Double, Point] with SimpleVectorStore[Double] {
//override newBuilder
}``````

Note that we also need a `VectorStore` to instantiate a vector. In this case, for simplicity we use a `SimpleVectorStore[Double]` which is the default implementation trait.

The method we have to override here, is `newBuilder` which we can implement from a static version to be defined later:

``override protected[this] def newBuilder: VectorBuilder[Double,Point] = Point.newBuilder``

Note the return type: `VectorBuilder[Double,Point]`, which is a vector builder that builds `Point`s from `Double` elements.

The remaining declarations go into the companion object:

``````object Point {
//define factory method
//define newBuilder
//define canBuildFrom
//define implicit class for scalar multiplication from left (optional)
}``````

First, since we made our constructor private, we need a factory method (as usual):

``def apply (x: Double, y: Double) = new Point(x, y)``

Then, we define the static `newBuilder` method mentioned before:

``````def newBuilder: VectorBuilder[Double,Point] =
new VectorBuilder[Double, Point] {
def result: Point = {
require(elems.length==2)
Point(elems(0), elems(1))
}
}``````

We use an anonymous implementation of the `VectorBuilder` trait. The trait provides an `elems` array which holds the elements put into the builder, so we just have to implement the `result` method to build a `Point` from the elements. Since our `Point` consists of `x` and `y` coordinates, we require that exactly 2 elements are there when building the result.

Third, we have to provide an implicit `VectorCanBuildFrom` implementation which follows the rules of the standard `scala.collection.generic.CanBuildFrom[From, Elem, To]` trait (i. e., define how to create a `To` collection from a `From` collection with elements of type `Elem`):

``````implicit def canBuildFrom: VectorCanBuildFrom[Point, Double, Point] =
new VectorCanBuildFrom[Point, Double, Point] {
def apply (): VectorBuilder[Double,Point] = newBuilder
def apply (from: Point): VectorBuilder[Double,Point] = newBuilder
}``````

Here, we just reuse our `newBuilder`.

Last, to return a `Point` for scalar multiplication from left (which is special as mentioned above), we define this simple implicit class:

``````implicit class ScalarOps (s: Double) {
def * (p: Point): Point = p * s
}``````

This last definition is optional, but nice.

The complete example is part of the `VectorSubtypeTest` class to be found under the `src/test/scala` folder. There is also a `MatrixSubtypeTest` with a similar example for subtyping the matrix trait.

### L6-4 Storage

Since version 2.0.0 vectors and matrices each use a separate storage engine that can be overridden. To this end, a vector inherits from a `VectorStore` definition trait which in turn has to be implemented by some implementation trait that realizes the storage; a similar mechanism works for matrices where the definition trait is called `MatrixStore`. The default implementation traits `SimpleVectorStore` (which uses `scala.collection.immutable.Vector` as storage for the vector elements) and `RowMatrixStore` (which uses a similar storage for the row vectors of a matrix) are mixed in during the standard creation process of a vector or matrix; see for instance the `DoubleVector` companion object to find out where this happens.

To demonstrate how this mechanism can be used to implement an own vector store for a subtype, we retouch the point example from the previous section. There we utilized the default `SimpleVectorStore[Double]` to store the elements of our point. In fact, this is a duplication of storage, since we already have the fields `x` and `y`.

So, here is an alternative that makes a point to be its own vector store:

``````class PointOSt private (val x: Double, val y: Double) extends
DoubleVector(1, Seq(x,y)) with IterableLike[Double, PointOSt] with
VectorLike[Double, PointOSt] {
//define index
//define apply
//define hash code
//override newBuilder
}
``````

Compared to the original `Point`, we left out the mix in of `SimpleVectorStore[Double]` and added three more definitions instead.

At first, a vector needs an `index` range. This is simple in this case, because we only have two elements numbered from 1 to 2:

``val index: Index = Index(1, 2)``

Second, we must define an `apply` method that maps an index to its value:

``````def apply (i: Int): Double = i match {
case 1 ⇒ x
case 2 ⇒ y
case _ ⇒ 0.0
}``````

Note that we mapped each integer other than 1 or 2 to zero, following the LinAlg concept of virtual elements.

The third definition we have to add is a `dataHashCode` that provides a hash code of the stored data elements (`x` and `y` in this case). Here, we don't need to take the index range into account, since this is done automatically when an overall vector hash code is computed. So, the hash-code implementation could, e. g., look like this:

``````protected val dataHashCode: Int = {
val prime = 31
var result = 1
result = prime*result+x.hashCode()
result = prime*result+y.hashCode()
result.toInt
}``````

We already had overridden the `newBuilder` method in the original `Point` example. We do not repeat the code here, since the only change is to replace the `Point` type by `PointOSt`; the same is true for the companion object that we can also copy from there with the change mentioned.

Beware that it is important to override `newBuilder` and to define an implicit `canBuildFrom` in the companion object, because otherwise the results of operations (like `p+q` for `p` and `q` of type `PointOSt`) would not use our own storage.

The complete example is part of the `VectorSubtypeOwnStoreTest` class to be found under the `src/test/scala` folder. There is also a `MatrixSubtypeOwnStoreTest` with a similar example for overriding the matrix storage engine.

## A Appendix

### Maven Coordinates

See The Central Repository. Choose the version you want (of course, the latest one is recommended) and look under "Dependency Information".

LinAlg - Scala Library for Vector and Matrix Types and Operations

Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

### Major Changes

Here we are describing compatibility-breaking changes.

#### Version 3.0.0

• New package structure: building, factory and storage classes and traits are moved to respective sub packages. Implies that import statements for that classes and traits have to be adapted.

• Zero vectors and matrices are now always equal regardless indices.

• Shorten on zero vectors and matrices now disregards indices and always yields lower index 1.

• The `VectorFactory` trait has an additional abstract `create` method.

#### Version 2.0.0

• Vector-builder class replaced by trait with target type parameter. In existing code, replace `new VectorBuilder...` by `VectorBuilder...` or even better `Vector.newBuilder...`.

• Matrix-builder class replaced by trait with source and target type parameter. In existing code, replace `new MatrixBuilder...` by `MatrixBuilder...` or even better `Matrix.newBuilder...`.

• Made `newbuilder` methods more specific to return a vector or matrix builder rather than a general builder.

• Specific `CanBuildFrom` traits introduced that return a vector or matrix builder rather than a general builder.

• Vector and matrix arithmetics moved to separate implementation traits with target type parameter. This allows subtypes of `Vector[E]` and `Matrix[E}` to get proper result types of operations by overriding newBuilder and implementing a suitable `CanBuildFrom` implicit.

• Storage-related stuff moved into separate traits, allowing subtypes to implement their own storage strategy. The implementation classes `Double/Int/NumericVector` and `RowMatrix` are abstract now. The default storage strategy is mixed in by the companion objects's apply method, respectively. As a consequence, the runtime class of a vector is not simply `Double/Int/NumericVector` anymore, but has the storage-implementation trait mixed in; likewise, the runtime class of a matrix is `RowMatrix` plus storage-implementation trait now.

### References

[1] S. G. van der Meulen, P. Kühling, "Programmieren in ALGOL68", Bd. II (Berlin, New York: de Gruyter), 149-188 (1977).