Feature request: symmetric matrix constructors

[ghost]

Some common and useful symmetric matrices it would be nice to construct in one step:

• Toeplitz
• Hermitian
• Circulant

https://en.wikipedia.org/wiki/Symmetric_matrix

[ghost]

Using the function I mentioned in #145, here's a sketch of a symmetric matrix constructor:

fn upper_diag_inds(nrows: usize, ncols: usize) -> Vec<[usize; 2]> {
  (0..nrows)
    .enumerate()
    .flat_map(|(iter_ind, row_ind)| {
      (iter_ind..ncols).map(move |col_ind| {
        [row_ind, col_ind]
      })
    }).collect()
}

// Computes the inverse of the triangular number function. Returns Result,
// because not every number is triangular.
fn triangular_inverse(tri_num: usize) -> Result<usize, Error> {
  let test_num = (((tri_num * 8 + 1) as Float).sqrt() - 1) / 2;
  if test_num == test_num.floor() {
    Ok(test_num)
  } else {
    Err("Input {} was not triangular.", &tri_num);
  }
}

// Assumes that the values are passed in row-major order, with all values below
// the diagonal omitted.
fn symmetric_matrix<T>(values: Vec<T>) -> Matrix<T>
  where T: NumCast + Float + Any {

  // Panics here if the number of values for the upper diagonal is not a
  // triangular number.
  let dim: usize = triangular_inverse(values.len()).unwrap();

  let mat: Matrix<T> = Matrix::zeros(dim, dim);

  for (val_ind, upper_diag_ind) in upper_diag_inds(dim, dim).enumerate() {
    mat[upper_diag_ind] = values[val_ind];

    // Add elements to transpose
    if upper_diag_ind[0] != upper_diag_ind[1] {
      mat[upper_diag_ind[upper_diag_ind[1],
                         upper_diag_ind[0]]] = values[val_ind];
    }
  }
  mat
}

[AtheMathmo]

Sorry for my slow response on this!

I agree with you that adding simple constructors for these matrices could be a good thing. I'll give my thoughts on this issue and #145 together here.

Firstly I think that because these are fairly niche things they are probably better kept a little apart from the current Matrix API. Additionally we have had some discussion about Symmetric/Triangular matrix storage in the past and these types seem like they may be better tackled as part of that work. See @Andlon 's work on the Permutation matrix for an idea of what I'm talking about.

As for the upper/lower_diag_inds my feeling is that these do not serve a common use case and it's a little difficult to understand exactly when they are useful. For example why not just do:

for i in 0..rows {
    for j in i..cols {

    }
}

As always I'm happy to hear thoughts from the other side. Let me know if I am missing anything!

[ghost]

I'll think more deeply about the other matrix items. I think you're right that there is potential if matrices that are known to have special properties are given "special" status within the library, certainly if they are known to be Symmetric.

The reason why it's useful to have a Vec or Iterator expression of the symmetric indices as opposed to just a couple of nested for loops is I think a standard characteristic of Rust: every iterative operation has an expressive functional equivalent. I don't think there is a one-liner that performs the task of creating an Iterator or Range like (0..rows).zip(cur_row..cols), and even if there was, it's more helpful for it to have a name.

On the other hand, if there was a SymmetricMatrix trait or struct, such an Iterator would follow naturally. What do you think?

[Andlon]

I think in general it's nice to provide some functionality for working with commonly occurring matrices with special structure. However, I think we should maybe think about this from a usage perspective. What applications actually require this functionality?

Let's take the symmetric matrix construction example. The proposed function takes a vector of values corresponding to the upper triangular part of the matrix. In most applications that I can think of, it's significantly harder to build a vector of such values than to work with a full matrix, and such a function would just copy the values into a full matrix anyway, so it would be both more intuitive and more efficient to just build a full matrix in the first place, using a double for loop or whatever is appropriate for your application.

A useful alternative, as @AtheMathmo mentioned, would be to have a custom struct for SymmetricMatrix which actually stores only the symmetric part of the matrix. This would save roughly half the memory costs of storing a full matrix. However, the development costs would be very high: we would have to provide specialized operator overloads (multiplication, add, subtract etc.) between all types of such special matrices for it to be useful. In some cases this might be worth it (in the case of PermutationMatrix, significant memory and performance benefits were attained in some cases), but the bar is quite high.

Now, the Toeplitz matrix is interesting, however. Like the PermutationMatrix, its special structure allows for much more efficient storage and algoritms. One can solve an invertible Toeplitz matrix in only O(n²) operations, and it can be stored in only O(n) memory. Moreover, we wouldn't really need to implement binary operators (I think? Does one often multiply by a Toeplitz matrix? I think not. And in any case it can be converted to Matrix for those purpose, although it is not optimal), so it wouldn't have to be very invasive to implement.

Finally, I just want to say that while I think custom storage for e.g. triangular, symmetric and other types of matrices is too much work to be worth it, a more enticing approach is the concept of special matrix "wrappers" through which one could access special APIs for special matrices. We've discussed this before in various different issues. As an example, one could imagine this (purely imaginary, not sure if it would work out this way):

let x = matrix![ 1.0, 0.0, 0.0;
                 2.0, 1.0, 0.0;
                 1.0, 2.0, 3.0];

// Wrapper verifies that the matrix satisfies the lower triangular constraint
let lower_triangular = LowerTriangular::from_matrix(x).expect("Matrix is lower triangular");
// Wrapper makes no verification
let lower_triangular = LowerTriangular::from_matrix_unchecked(x);
// Wrapper verifies using an element-wise tolerance rather than exactly zero
let lower_triangular = LowerTriangular::from_matrix_with_tol(x, 1e-12)
                                      .expect("Matrix is lower triangular to given precision");

// Once obtained, one can efficiently solve the triangular system
let y = lower_triangular.solve(Vector::zeros(3))
                                    .expect("Triangular system is invertible.")

// Or compute the determinant
let det = lower_triangular.det();

A major benefit of such an approach, is that it allows us to check that the constraint is satisfied exactly once, and then solve many such systems without risk of the constraint breaking between different calls to solve.