Idea: Add immutable diagonal matrix view with O(n) space

[Waterdragen]

Current state

Currently, ndarray provides:
Array2::from_diag: Create a mutable 2D array from a vector
Array2::from_diag_elem: Create a mutable 2D array with repeated diagonal elements
Array.diag(): Make a view of the diagonal elements of the ndarray

Problem

However, the mutability might not be needed, in this case, only the diagonal elements are non-zero, thus we only need O(n) space. Using a large Array might waste a lot of space.

Possible use case

Let's call the new struct DiagArray, and since it can represent diagonal arrays of any dimension, we have to specify the dimensions. Here is a possible use case:

let my_array = Array1::from_vec(vec![2, 3, 4]);
let weights = DiagArray::from_vec(2, vec![2, 3, 4]);

assert_eq!(weights[[0, 0]], 2);

let my_array_squared = my_array.dot(&weights);
assert_eq!(my_array_squared, Array1::from_vec(vec![4, 9, 16]));

Unresolved questions

• Since only Array2 has from_diag, do we really need diagonal arrays for all dimensions?
• Should we use type name or integer parameter to specify the dimensions?

[akern40]

Thank you for the idea! I am particularly curious to know whether you ran across a use case for this in your own work, or whether this is driven more by curiosity.

Right now, ndarray is really engineered for dense arrays and centered around just a few types. There's been discussion - and I have some ongoing work - that would move us towards a more general, trait-based approach. But for the moment, implementing this kind of sparse diagonal matrix would be very difficult.

[Waterdragen]

Hi there! I use diagonal matrices mainly for weighted least squares adjustment. The formula looks like this:

$\hat{\mathit{\beta }}=({\mathbf{A}}^{\mathrm{\top }}\cdot \mathbf{P}\cdot \mathbf{A}{)}^{-1}\cdot {\mathbf{A}}^{\mathrm{\top }}\cdot \mathbf{P}\cdot \mathbf{l}$

The weight matrix P is a diagonal matrix. Here is a real-world example:

and my solution to question 1 adjusted height

use ndarray::linalg::Dot;
use ndarray::{Array1, Array2};
use ndarray_linalg::solve::Inverse;
use sprs::{CsMat, TriMat};

fn weighted_least_squares(a: &Array2<f64>, p: &CsMat<f64>, l: &Array1<f64>) -> Array1<f64> {
    a.t().dot(p).dot(a).inv().unwrap().dot(&a.t()).dot(p).dot(l)
}

fn main() {
    let a: Array2<f64> = Array2::from_shape_vec((6, 3), vec![
     //  B    C    D
         1.,  0.,  0.,
        -1.,  1.,  0.,
         0., -1.,  1.,
         0.,  0., -1.,
        -1.,  0.,  1.,
         0.,  1.,  0.
    ]).unwrap();

    let mut p = TriMat::new((6, 6)); // 6x6 triangular matrix to be converted to CsMat
    p.add_triplet(0, 0, 1./36.); // Diagonal element at (0, 0)
    p.add_triplet(1, 1, 1./16.); // Diagonal element at (1, 1)
    p.add_triplet(2, 2, 1./25.); // Diagonal element at (2, 2)
    p.add_triplet(3, 3, 1./9.);  // Diagonal element at (3, 3)
    p.add_triplet(4, 4, 1./16.); // Diagonal element at (4, 4)
    p.add_triplet(5, 5, 1./144.);  // Diagonal element at (5, 5)
    let p: CsMat<f64> = p.to_csr();
    // OR with `ndarray::Array2`
    // let p_alt = Array2::from_diag(&Array1::from_vec(vec![
    //     1./36., 1./16., 1./25., 1./9., 1./16., 1./144.
    // ]));

    let l: Array1<f64> = Array1::from_vec(vec![10.509, 5.360, -8.523, -7.348, -3.167, 15.881]);

    let mut heights = weighted_least_squares(&a, &p, &l);

    let height_a= 437.596;
    heights += height_a;
    let (height_b, height_c, height_d) = (heights[0], heights[1], heights[2]);
    assert_eq!(height_a, 437.596);
    assert_eq!(height_b, 448.1087117287866);
    assert_eq!(height_c, 453.4684677834756);
    assert_eq!(height_d, 444.94360533133613);

    println!("{heights}");
}