Linear Algebra

Adding matrices

Creates two 2-D matrices with ndarray::arr2 and sums them element-wise.

Note the sum is computed as let sum = &a + &b. The & operator is used to avoid consuming a and b, making them available later for display. A new array is created containing their sum.

use ndarray::arr2;

fn main() {
    let a = arr2(&[[1, 2, 3],
                   [4, 5, 6]]);

    let b = arr2(&[[6, 5, 4],
                   [3, 2, 1]]);

    let sum = &a + &b;

    println!("{}", a);
    println!("+");
    println!("{}", b);
    println!("=");
    println!("{}", sum);
}

Multiplying matrices

Creates two matrices with ndarray::arr2 and performs matrix multiplication on them with ndarray::ArrayBase::dot.

use ndarray::arr2;

fn main() {
    let a = arr2(&[[1, 2, 3],
                   [4, 5, 6]]);

    let b = arr2(&[[6, 3],
                   [5, 2],
                   [4, 1]]);

    println!("{}", a.dot(&b));
}

Multiply a scalar with a vector with a matrix

Creates a 1-D array (vector) with ndarray::arr1 and a 2-D array (matrix) with ndarray::arr2.

First, a scalar is multiplied by the vector to get another vector. Then, the matrix is multiplied by the new vector with ndarray::Array2::dot. (Matrix multiplication is performed using dot, while the * operator performs element-wise multiplication.)

In ndarray, 1-D arrays can be interpreted as either row or column vectors depending on context. If representing the orientation of a vector is important, a 2-D array with one row or one column must be used instead. In this example, the vector is a 1-D array on the right-hand side, so dot handles it as a column vector.

use ndarray::{arr1, arr2, Array1};

fn main() {
    let scalar = 4;

    let vector = arr1(&[1, 2, 3]);

    let matrix = arr2(&[[4, 5, 6],
                        [7, 8, 9]]);

    let new_vector: Array1<_> = scalar * vector;
    println!("{}", new_vector);

    let new_matrix = matrix.dot(&new_vector);
    println!("{}", new_matrix);
}

Vector comparison

The ndarray crate supports a number of ways to create arrays -- this recipe creates ndarray::Arrays from std::Vec using from. Then, it sums the arrays element-wise.

This recipe contains an example of comparing two floating-point vectors element-wise. Floating-point numbers are often stored inexactly, making exact comparisons difficult. However, the assert_abs_diff_eq! macro from the approx crate allows for convenient element-wise comparisons. To use the approx crate with ndarray, the approx feature must be added to the ndarray dependency in Cargo.toml. For example, ndarray = { version = "0.13", features = ["approx"] }.

This recipe also contains additional ownership examples. Here, let z = a + b consumes a and b, updates a with the result, then moves ownership to z. Alternatively, let w = &c + &d creates a new vector without consuming c or d, allowing their modification later. See Binary Operators With Two Arrays for additional detail.

use approx::assert_abs_diff_eq;
use ndarray::Array;

fn main() {
  let a = Array::from(vec![1., 2., 3., 4., 5.]);
  let b = Array::from(vec![5., 4., 3., 2., 1.]);
  let mut c = Array::from(vec![1., 2., 3., 4., 5.]);
  let mut d = Array::from(vec![5., 4., 3., 2., 1.]);

  let z = a + b;
  let w =  &c + &d;

  assert_abs_diff_eq!(z, Array::from(vec![6., 6., 6., 6., 6.]));

  println!("c = {}", c);
  c[0] = 10.;
  d[1] = 10.;

  assert_abs_diff_eq!(w, Array::from(vec![6., 6., 6., 6., 6.]));

}

Vector norm

This recipe demonstrates use of the Array1 type, ArrayView1 type, fold method, and dot method in computing the l1 and l2 norms of a given vector. + The l2_norm function is the simpler of the two, as it computes the square root of the dot product of a vector with itself. + The l1_norm function is computed by a fold operation that sums the absolute values of the elements. (This could also be performed with x.mapv(f64::abs).scalar_sum(), but that would allocate a new array for the result of the mapv.)

Note that both l1_norm and l2_norm take the ArrayView1 type. This recipe considers vector norms, so the norm functions only need to accept one-dimensional views (hence ArrayView1). While the functions could take a parameter of type &Array1<f64> instead, that would require the caller to have a reference to an owned array, which is more restrictive than just having access to a view (since a view can be created from any array or view, not just an owned array).

Array and ArrayView are both type aliases for ArrayBase. So, the most general argument type for the caller would be &ArrayBase<S, Ix1> where S: Data, because then the caller could use &array or &view instead of x.view(). If the function is part of a public API, that may be a better choice for the benefit of users. For internal functions, the more concise ArrayView1<f64> may be preferable.

use ndarray::{array, Array1, ArrayView1};

fn l1_norm(x: ArrayView1<f64>) -> f64 {
    x.fold(0., |acc, elem| acc + elem.abs())
}

fn l2_norm(x: ArrayView1<f64>) -> f64 {
    x.dot(&x).sqrt()
}

fn normalize(mut x: Array1<f64>) -> Array1<f64> {
    let norm = l2_norm(x.view());
    x.mapv_inplace(|e| e/norm);
    x
}

fn main() {
    let x = array![1., 2., 3., 4., 5.];
    println!("||x||_2 = {}", l2_norm(x.view()));
    println!("||x||_1 = {}", l1_norm(x.view()));
    println!("Normalizing x yields {:?}", normalize(x));
}

Invert matrix

Creates a 3x3 matrix with nalgebra::Matrix3 and inverts it, if possible.

use nalgebra::Matrix3;

fn main() {
    let m1 = Matrix3::new(2.0, 1.0, 1.0, 3.0, 2.0, 1.0, 2.0, 1.0, 2.0);
    println!("m1 = {}", m1);
    match m1.try_inverse() {
        Some(inv) => {
            println!("The inverse of m1 is: {}", inv);
        }
        None => {
            println!("m1 is not invertible!");
        }
    }
}

(De)-Serialize a Matrix

Serialize and deserialize a matrix to and from JSON. Serialization is taken care of by serde_json::to_string and serde_json::from_str performs deserialization.

Note that serialization followed by deserialization gives back the original matrix.

use nalgebra::DMatrix;

fn main() -> Result<(), std::io::Error> {
    let row_slice: Vec<i32> = (1..5001).collect();
    let matrix = DMatrix::from_row_slice(50, 100, &row_slice);

    // serialize matrix
    let serialized_matrix = serde_json::to_string(&matrix)?;

    // deserialize matrix
    let deserialized_matrix: DMatrix<i32> = serde_json::from_str(&serialized_matrix)?;

    // verify that `deserialized_matrix` is equal to `matrix`
    assert!(deserialized_matrix == matrix);

    Ok(())
}