diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs
index 182b8a6d..bce2ea41 100644
--- a/src/linalg/udu.rs
+++ b/src/linalg/udu.rs
@@ -4,10 +4,21 @@ use serde::{Deserialize, Serialize};
 use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, MatrixN, VectorN, U1};
 use crate::dimension::Dim;
+use crate::storage::Storage;
 use simba::scalar::ComplexField;
 
 /// UDU factorization
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(bound(serialize = "VectorN<N, D>: Serialize, MatrixN<N, D>: Serialize"))
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(bound(
+        deserialize = "VectorN<N, D>: Deserialize<'de>, MatrixN<N, D>: Deserialize<'de>"
+    ))
+)]
 #[derive(Clone, Debug)]
 pub struct UDU<N: ComplexField, D: Dim>
 where
@@ -36,33 +47,30 @@ where
     /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360
     pub fn new(p: MatrixN<N, D>) -> Self {
         let n = p.ncols();
-        let n_as_dim = D::from_usize(n);
+        let n_dim = p.data.shape().1;
 
-        let mut d = VectorN::<N, D>::zeros_generic(n_as_dim, U1);
-        let mut u = MatrixN::<N, D>::zeros_generic(n_as_dim, n_as_dim);
+        let mut d = VectorN::zeros_generic(n_dim, U1);
+        let mut u = MatrixN::zeros_generic(n_dim, n_dim);
 
         d[n - 1] = p[(n - 1, n - 1)];
-        u[(n - 1, n - 1)] = N::one();
-
-        for j in (0..n - 1).rev() {
-            u[(j, n - 1)] = p[(j, n - 1)] / d[n - 1];
-        }
+        u.column_mut(n - 1)
+            .axpy(N::one() / d[n - 1], &p.column(n - 1), N::zero());
 
         for j in (0..n - 1).rev() {
+            let mut d_j = d[j];
             for k in j + 1..n {
-                d[j] = d[j] + d[k] * u[(j, k)].powi(2);
+                d_j += d[k] * u[(j, k)].powi(2);
             }
 
-            d[j] = p[(j, j)] - d[j];
+            d[j] = p[(j, j)] - d_j;
 
             for i in (0..=j).rev() {
+                let mut u_ij = u[(i, j)];
                 for k in j + 1..n {
-                    u[(i, j)] = u[(i, j)] + d[k] * u[(j, k)] * u[(i, k)];
+                    u_ij += d[k] * u[(j, k)] * u[(i, k)];
                 }
 
-                u[(i, j)] = p[(i, j)] - u[(i, j)];
-
-                u[(i, j)] /= d[j];
+                u[(i, j)] = (p[(i, j)] - u_ij) / d[j];
             }
 
             u[(j, j)] = N::one();