get jacobian and hessian not only row by row

[eweca-d]

I am using this crate to do curve fitting with 'nalgebra' that is an algebra crate. For all I know, most algebra library tends to use column-major matrix and store matrix in a vec storage. Hence, I think it is better to provide new api or pass the paramter to control the output matrix shape to be:

• row-major: [[T; ncols]; nrows]
• column-major: [[T; nrows]; ncols]
• contiguous column by column: [T; nrows * ncols] => default for many crates
• contiguous row by row: [T; nrows * ncols]

Thank you for your great job. It works very well in my crate right now.

[kmolan]

Hi @eweca-d I see only three possible output matrices. What's the difference between contiguous column by column and contiguous row by row?

I will need some time to think of the implementation. What you're suggesting is helpful, but would require separate functions since the output type is for each case is different. Ideally I would like one function that returns the output no matter what the desired output shape is. Doing such things with vectors is easy, but also uses heap allocations which I want to avoid.

Give me some to think about it. I will update you soon.

[eweca-d]

"contiguous column by column" is the flatten array of column-major matrix, "contiguous row by row" is the flatten array of row-major matrix. There can be constructed by:

pub fn get_custom<T: ComplexFloat, const NUM_FUNCS: usize, const NUM_VARS: usize>(function_matrix: &[&dyn Fn(&[T; NUM_VARS]) -> T; NUM_FUNCS], vector_of_points: &[T; NUM_VARS], step_size: f64, mode: mode::DiffMode) -> Result<[[T; NUM_VARS]; NUM_FUNCS], ErrorCode>
{
    if function_matrix.is_empty()
    {
        return Err(ErrorCode::VectorOfFunctionsCannotBeEmpty);
    }

    // row-major
    let mut result = [[T::zero(); NUM_VARS]; NUM_FUNCS];
    for row_index in 0..NUM_FUNCS
    {
        for col_index in 0..NUM_VARS
        {
            result[row_index][col_index] = single_derivative::get_partial_custom(&function_matrix[row_index], col_index, vector_of_points, step_size, mode)?;
        }
    }

    // col-major
    let mut result = [[T::zero(); NUM_FUNCS]; NUM_VARS];
    for col_index in 0..NUM_VARS
    {
        for row_index in 0..NUM_FUNCS
        {
            result[col_index][row_index] = single_derivative::get_partial_custom(&function_matrix[row_index], col_index, vector_of_points, step_size, mode)?;
        }
    }

    // contiguous row by row
    let mut result = [T::zero(); NUM_VARS * NUM_FUNCS];
    for row_index in 0..NUM_FUNCS
    {
        for col_index in 0..NUM_VARS
        {
            result[row_index * NUM_VARS + col_index] = single_derivative::get_partial_custom(&function_matrix[row_index], col_index, vector_of_points, step_size, mode)?;
        }
    }

    // contiguous col by col
    let mut result = [T::zero(); NUM_FUNCS * NUM_VARS];
    for col_index in 0..NUM_VARS
    {
        for row_index in 0..NUM_FUNCS
        {
            result[col_index * NUM_FUNCS + row_index] = single_derivative::get_partial_custom(&function_matrix[row_index], col_index, vector_of_points, step_size, mode)?;
        }
    }
    
    return Ok(result);
}

---

As for the heap allocations, I do not think it is a problem. The current implement can just work for a small size Jacobian matrix. However, if I try to fit a curve with 20000 data points, the stack will overflow. Furthermore, the vector has no significant impact on the performance if the size is relatively large. I recommend replace the array of array by nalgebra Matrix. It provides static allocated and dynamic matrix for the use of both small and large matrix.

[eweca-d]

We can return [T; N * M] in the custom function. Then, we can convert it to the shape we desired without reallocation.

An example:

fn main() {    
    let a = [0.0, 0.1, 0.2];
    let b = [1.1, 1.2];
    let jacobian = get_jacobian::<3, 2, 6>(a, b);
    println!("{:?}", jacobian);
}

fn get_jacobian_custom<const N: usize, const M: usize, const L: usize>(array1: [f64; N], array2: [f64; M]) -> [f64; L] {
    let mut jacobian = [0.0; L];
    for i in 0..N {
        for j in 0..M {
            jacobian[i * M + j] = array1[i] * array2[j];
        }
    }
    jacobian
}

fn get_jacobian<const N: usize, const M: usize, const L: usize>(array1: [f64; N], array2: [f64; M]) -> [[f64; M]; N] {
    let jacobian: [f64; L] = get_jacobian_custom(array1, array2);
    let ptr = jacobian.as_ptr() as *const [[f64; M]; N];
    unsafe { *ptr }
}

We can only replace L with N * M directly in nightly channel, but I guess we could find a workaround to solve it to eliminate the annoyed generic annotation.

[kmolan]

You point about stack overflow is a good one. Unfortunately even with nalgebra the static matrices (SMatrix) would soon overflow with 20000 points. What we really need here is a DMatrix for heap allocation.

I have started a bigger task of providing methods using vectors and DMatrices as an optional feature of the crate. That one might take some time, but I'm working on it!

[eweca-d]

Thank you for your effort! I really expect to enable the new features in my work.

[kmolan]

https://github.com/kmolan/multicalc-rust?tab=readme-ov-file#16-experimental
Chose to keep the matrix shape the same for API compatibility. Tried nalgebra DMatrix but num-complex traits don't play nice with it (they don't implement the Debug trait).

Long term solution is to hand-roll my own numeric trait, but for now this should atleast help you with stack overflows.

[eweca-d]

Thank you so much! It will really help if I derive a big Jacobian matrix in my work!