SVD on SMatrix with const generic size #1480
Replies: 1 comment
-
|
Hey, I don't have much experience with generic expressions and as they are experimental, I chose a different option. I hope it is still interesting for your use-case :) I mainly follow the approach from nalgebra's documentation. Unfortunately, the website is not up-to-date but you find a newer version here. We use the feature build around the Dim trait from nalgebra. Dim is implemented by Const and by Dyn in nalgebra. You can create a Matrix with Const and Dyn using OMatrix:
Neat right? This allows us to write code that either uses Dyn or Const so we can decide later if we want static or dynamic objects / stack or heap allocation, e.g.: use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, OMatrix, Scalar, Dyn, Const};
fn generic_function<T: Scalar, R: Dim, C: Dim>(input: OMatrix<T, R, C>)
where
DefaultAllocator: nalgebra::allocator::Allocator<R, C> // Important! I believe we need this so we are able to allocate the matrix.
{
todo!()
}
fn main() {
let x = OMatrix::<f64, Const<3>, Dyn>::zeros(4);
generic_function(x);
}Basically the major downside of the approach compared to the const generic expression is that we have to carry a lot of traits with us. Thiy is anoying but it is not too difficult :) This same holds when using operators such as +, *, ... etc. Instead of using the opperators directly nalegbra includes DimSum, DimProd, ... So instead of writing fn main() {
let x = OMatrix::<f64, Const<3>, Dyn>::zeros(4);
let (n_rows, n_columns) = x.shape_generic();
let extended_matrix = OMatrix::<f64, _, _>::zeros_generic(DimAdd::add(n_rows, Const::<1>), n_columns);
}Since nalegbra uses pretty cool but hacky trait stuff, the story is a bit more complicated if we look at generics. If we have a generic use nalgebra::{DefaultAllocator, DimAdd, DimSum, allocator::Allocator, Dim, OMatrix, Scalar, Dyn, Const};
fn extend_columns_by_one<T: Scalar, R: Dim, C: Dim>(input: OMatrix<T, R, C>)
where
C: DimAdd<Const<1>>,
DefaultAllocator:
nalgebra::allocator::Allocator<R, C>
+ nalgebra::allocator::Allocator<R, DimSum<C, Const<1>>>
{
let (n_rows, n_columns) = input.shape_generic();
let extended_matrix = OMatrix::<f64, _, _>::zeros_generic(n_rows, DimAdd::add(n_columns, Const::<1>));
}
fn main() {
let x = OMatrix::<f64, Const<3>, Dyn>::zeros(4);
extend_columns_by_one(x);
}Please take this with a grain of salt, since I figured that stuff out on my own, so probably not all is 100% correct. I find it so cool that they basically used the trait system to implement const generic impressions by themself :) Putting all together gives us: use nalgebra::{allocator::Allocator, Const, DefaultAllocator, Dim, DimDiff, DimMin, DimMinimum, DimMul, DimName, DimProd, DimSub, OMatrix, OVector, SMatrix, SVector, Scalar, ToTypenum, U1, U11, U8};
use pyo3_stub_gen::derive;
trait CEquationparametersAllocator<D: Dim>: Allocator<D, Const<6>> + Allocator<D, Const<2>> + Allocator<D, Const<4>> {}
struct CEquationparameters<T: Scalar, D: Dim> where
DefaultAllocator: CEquationparametersAllocator<D>
{
points: OMatrix<T, D, Const<6>>,
weights: OMatrix<T, D, Const<2>>,
coefficients: OMatrix<T, D, Const<4>>
}
// Okay lets look at the trait bounds of the svd method:
// impl<T, R, C, S> Matrix<T, R, C, S>
// pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<T, R, C>
// where
// R: DimMin<C>,
// DimMinimum<R, C>: DimSub<U1>,
// DefaultAllocator: Allocator<R, C> + Allocator<C> + Allocator<R>
// + Allocator<DimDiff<DimMinimum<R, C>, U1>> + Allocator<DimMinimum<R, C>, C> + Allocator<R, DimMinimum<R, C>> + Allocator<DimMinimum<R, C>>,
// ...
// Thats alot of Allocators! >_< Lets bundle them to make it more readable :')
// Lets bundle the first line for the DefaultAllocators:
trait SVDAllocator<R: Dim, C: Dim>: Allocator<R, C> + Allocator<C> + Allocator<R>
{}
// And lets bundle the second line for the DefaultAllocator:
trait SVDDIMMinimumAllocator<D: Dim, R: Dim, C: Dim> : Allocator<DimDiff<D, U1>> + Allocator<D, C> + Allocator<R, D> + Allocator<D>
where
D: DimSub<U1>
{}
impl<T: Scalar, D: DimName> CEquationparameters<T, D> where
// we need N*11 and N*8 to work, thats done by DimMul
D: DimMul<U11> + DimMul<U8>, // now we can multiply instead of writing N*8 we use DimProd<D, U8>
// Further we need to be able to create/allocate a OVector<N*11>
DefaultAllocator: Allocator<DimProd<D, U11> >, // allocator for x_0 vector
// An we need an allcoator for the Matrix<N*8, N*11>
DefaultAllocator: Allocator<DimProd<D,U8>, DimProd<D,U11>>, // allocator for derivative matrix
// We introduce the allocator for the CEquationParameter
DefaultAllocator: CEquationparametersAllocator<D>,
// Finally we have to introduce all the SVD traits, since we apply SVD on the derivative matrix we have:
// R = DimProd<D,U8>
// C = DimProd<D,U11>
// Okay lets get this over with, lets copy past this from the SVD deffinition above and replace R and C
DimProd<D,U8>: DimMin<DimProd<D,U11>>,
DimProd<D,U8>: DimMin<DimProd<D,U11>>,
DimMinimum<DimProd<D,U8>, DimProd<D,U11>>: DimSub<U1>, // see SVD trait bounds
DefaultAllocator: SVDAllocator<DimProd<D,U8>, DimProd<D,U11>>
+ SVDDIMMinimumAllocator<DimMinimum<DimProd<D,U8>, DimProd<D,U11>>,DimProd<D,U8>, DimProd<D,U11>>
{
fn to_vec(&self) -> OVector<f32, DimProd<D, U11> >
{
todo!()
}
fn derivative(&self) -> OMatrix<f32, DimProd<D, U8>, DimProd<D, U11>>
{
todo!()
}
pub fn improve_values(&self)
{
let derivative = self.derivative();
let x_0 = self.to_vec();
let b = (&derivative)*x_0;
let result = derivative.svd(true, true).solve(&b, 1e-6);
}
} |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
-
I have the following situation:
So I need to do a SVD on a statically sized matrix depending on a
const Nsize parameter, e.g. number of points to minimize an equation. With the above code, rustc rises the following error:Introducing the suggested error bound, so
will result in the following error:
However, introducing this bound will result in a cycle:
Is there any way to get the compile-time safety of the matrix calculation?
Beta Was this translation helpful? Give feedback.
All reactions