Skip to content
Cédric Boutillier edited this page Mar 19, 2014 · 10 revisions

NMatrix is a linear algebra library for Ruby, written mostly in C. It can be used with or without SciRuby, but is part of the SciRuby project. NMatrix was inspired by and based heavily upon NArray, by Masahiro Tanaka.

Installation

See Installation for details.

Generally, you'd only need gem install nmatrix. If you want to install from source:

git clone git://github.com/SciRuby/nmatrix.git
cd nmatrix
rake install

Features

In addition to dense matrices, NMatrix supports two types of sparse storage (list-of-list and Yale). We support the following data classes:

  • Byte (unsigned 8-bit integer)
  • Integer (8, 16, 32, and 64-bit)
  • Floating point (32 and 64-bit)
  • Complex (64 and 128-bit)
  • Rational (32, 64, and 128-bit)
  • Ruby objects

NMatrix makes use of and exposes the CBLAS API included in ATLAS. However, we have also written BLAS templates for NMatrix-supported data types that are not included in ATLAS or LAPACK (e.g., rationals).

Differences from other matrix libraries

In contrast to ruby's Matrix class, and scipy.matrix, the * operator is used for element-wise multiplication and the #dot method for matrix-matrix or matrix-vector multiplication.

Brief Tutorial

Creating matrices

Currently, you create NMatrix objects with the NMatrix#new method. The first parameter is the dimensions of the matrix, and the second parameter contains the initial values. For example, the following creates a 3x4 matrix filled with zeros:

>> n = NMatrix.new( [3, 4], 0)
0  0  0  0
0  0  0  0
0  0  0  0
=> nil

If the first parameter is an integer, NMatrix.new returns a square matrix.

>> n = NMatrix.new( 4, 0)   # Same as NMatrix.new( [4,4], 0 )
0  0  0  0
0  0  0  0
0  0  0  0
0  0  0  0
=> nil

Matrices need not have only two dimensions:

# Create a 3D matrix (4 x 4 x 4)
b = NMatrix.new([4,4,4], 0)

In addition, the second parameter to NMatrix.new can be an array which specifies all the values in the matrix:

>> B = NMatrix.new([2,5], [1,2,3,4,5,6,7,8,9,0])
1 2 3 4 5
6 7 8 9 0

If the array of values has fewer elements than the matrix requires, NMatrix.new will repeat the array as many times as required:

>> q = NMatrix.new([3,6], [1,2,3])
1  2  3  1  2  3
1  2  3  1  2  3
1  2  3  1  2  3

Simpler constructor for matrices

Thanks to Daniel Carrera, we have an in-line constructor that makes creating matrices a very simple and readable operation:

a = N[ [1,2,3,4] ]                    # =>  [ [1.0  2.0  3.0  4.0] ]
a = N[ [1,2,3,4], dtype: :int32 ]     # =>  [ [1    2    3    4] ]
a = N[ [1,2,3], [3,4,5] ]             # =>  1.0  2.0  3.0
                                      #     3.0  4.0  5.0

Note that N[1,2,3,4] will create a "dimensionless" (technically 1-dimensional) vector. We advise creating a vector with the same dimensionality as the matrices you're working with. For two dimensions:

N[ [1,2,3,4] ] # row vector
N[ [1],[2],[3],[4] ] # column vector
N[ [1,2,3,4] ].transpose # column vector

One-dimensional matrices

One-dimensional matrices aren't really one-dimensional. They are n-dimensional, but n-1 of their dimensions are length 1.

Here's how you create vectors with different orientations:

>> x = NMatrix.new([4], [1,2,3,4]) # no orientation (dimensionless/1-dimensional)
>> x = NMatrix.new([4,1], [1,2,3,4]) # column vector in two dimensions (4 rows, 1 column)
>> x = NMatrix.new([1,4], [1,2,3,4]) # row vector in two dimensions (1 row, 4 columns)
>> x = NMatrix.new([4,1,1], [1,2,3,4]) # column vector in three dimensions (4 rows, 1 column, 1 layer)

Matrix functions and operations

NMatrix supports the standard arithmetic operations on matrices:

# a = NMatrix.new([3,3], 2) # 3x3 matrix containing all 2's

a + a       # Element-wise / Matrix addition.
a - a       # Element-wise / Matrix subtraction.
a ** 3      # Element-wise exponentiation.
a / 2       # Element-wise division.
a * a       # Element-wise multiplication.
a.dot a     # Matrix dot-product multiplication.

Shortcuts

There are several shortcuts for creating NMatrix or NVector objects. They are as follows:

NMatrix.ones [3, 3]                     # A 3x3 matrix of ones.
NVector.ones 3                          # A 3-vector of ones.
NMatrix.zeros [4, 4] (or zeroes [4, 4]) # Creates a matrix of zeros with dimensions as parameters.
NVector.zeros 4                         # Creates a vector of zeros
NMatrix.identity [5, 5] (or eye [5, 5]) # A 5x5 identity matrix. Only works with NMatrix.
NMatrix.seq [3, 4]                      # A 3x4 matrix with values from 0 to 11.
NVector.seq 8                           # An 8-vector with a sequence of values from 0 to 7.
random [6,6]                            # A 6x6 matrix of random values between 0 and 1.

m.row(0)                                # Return the 1st column of the matrix m. Only works with NMatrix.
m.column(2)                             # Return the 3rd column of the matrix m. Only works with NMatrix.

NVector.indgen(8)                       # To make IDL users happy. Same as `NVector.seq(8, :int32)`.
NVector.findgen(8)                      # Same as `NVector.seq(8, :float32)`.
NVector.bindgen(8)                      # Same as `NVector.seq(8, :byte)`.
NVector.cindgen(8)                      # Same as `NVector.seq(8, :complex64)`.

NVector.linspace(0, pi, 100)            # A vector of 100 values from 0 to pi, inclusive. Only works with NVector.

Solving linear systems

For now, the best way is probably to use NMatrix::LAPACK::clapack_gesv, which is the Ruby-exposed API for ATLAS's CLAPACK implementation of SGESV/DGESV/CGESV/ZGESV.

FIXME: Needs documentation (and a spec)

Matrix Decomposition

Here is LU decomposition (using CLAPACK):

n = NMatrix.new(3, [4,9,2,3,5,7,8,1,6], dtype: :float64)
lu = n.factorize_lu # the result are the Lower and Upper components stored in a single matrix

There are multiple methods for singular value decomposition, all of which use LAPACK instead of CLAPACK -- so they should work even if you haven't installed your own ATLAS (e.g., on Mac OS X):

u, sigma, v_transpose = n.gesvd # LAPACK's GESVD function
u, sigma, v_transpose = n.gesdd # Divide-and-conquer strategy

You can also do eigenvalue decomposition using NMatrix::LAPACK::lapack_geev. See the specs for an example.

Some of these functions still need NMatrix instance methods. Please consider contributing one!

Advanced features

NMatrix supports several data types. NMatrix will try to guess the dtype based on the first initial value you provide (e.g., NMatrix.new(4, [0.0, 1]) will be :float32), but you can choose to provide a dtype in addition to or in lieu of initial values:

>> n = NMatrix.new(4, [0,1], dtype: :rational128)
0/1  1/1  0/1  1/1
0/1  1/1  0/1  1/1
0/1  1/1  0/1  1/1
0/1  1/1  0/1  1/1
>> m = NMatrix.new(4, dtype: :int64)  # Caution! No initialization of contents. Example output:
              8,  140642513449760,  0,                0
              0,                0,  0,                0
              0,                0,  0,  140642505889952
140642512627712,      -4294967290,  1,  140734747434816
>> m = NMatrix.new(4, dtype: :rational128) # not recommended (but allowed)
>> m = NMatrix.new(4) # 4x4 Ruby :object matrix containing only nil

The storage type (stype) can also be specified, prior to the dimension argument. However, with sparse storage formats, initial values don’t make sense, and these matrices will contain zeros by default.

# empty list-of-lists-of-lists 4x3x4 matrix
n = NMatrix.new([4,3,4], dtype: :int64, stype: :list)

# Ruby objects in a 'Yale' sparse matrix
m = NMatrix.new([5,4], stype: :yale)

# A byte matrix containing a gradient
o = NMatrix.new(5, [0,1,2,3,4], dtype: :byte)

The matrix m created above is a Yale-format sparse matrix, or more specifically, “new Yale,” which differs from “old Yale” in that the diagonal is stored separately from the non-diagonal elements. Thus, diagonals can be accessed and set in constant time.

Currently, all storage is row-based.

You can also convert between any of these three stypes using cast:

n = NMatrix.new(4, dtype: :int64, stype: :list)
n[0,0] = 5
n[0,3] = -2
dense = n.cast(:dense) # you may also supply other information to cast, e.g.,
dense_double = n.cast(stype: :dense, dtype: :float64)

Let's say you have a dense matrix consisting of only 1s with a few 0s. There are two ways to cast this to sparse matrices. You can have the sparse storage store only the non-zero values, or only the non-one values, by specifying a default value for cast.

list = dense.cast(stype: :list, default: 1) # 1-based matrix
list = dense.cast(:list) # default is 0-based

The same works for Yale matrices.

Clone this wiki locally