diff --git a/algorithms/linfa-clustering/Cargo.toml b/algorithms/linfa-clustering/Cargo.toml
index 4b034789f..a8340e159 100644
--- a/algorithms/linfa-clustering/Cargo.toml
+++ b/algorithms/linfa-clustering/Cargo.toml
@@ -70,3 +70,7 @@ harness = false
 [[bench]]
 name = "gaussian_mixture"
 harness = false
+
+[[bench]]
+name = "k_modes"
+harness = false
diff --git a/algorithms/linfa-clustering/benches/k_modes.rs b/algorithms/linfa-clustering/benches/k_modes.rs
new file mode 100644
index 000000000..e47f91e6d
--- /dev/null
+++ b/algorithms/linfa-clustering/benches/k_modes.rs
@@ -0,0 +1,104 @@
+use criterion::{
+    black_box, criterion_group, criterion_main, AxisScale, BenchmarkId, Criterion,
+    PlotConfiguration,
+};
+use linfa::benchmarks::config;
+use linfa::prelude::Fit;
+use linfa::DatasetBase;
+use linfa_clustering::{KModes, KModesInit};
+use ndarray::Array2;
+use ndarray_rand::rand::Rng;
+use ndarray_rand::rand::SeedableRng;
+use rand_xoshiro::Xoshiro256Plus;
+
+fn generate_synthetic_categorical(
+    n_samples: usize,
+    n_features: usize,
+    n_categories: usize,
+    rng: &mut Xoshiro256Plus,
+) -> Array2<usize> {
+    Array2::from_shape_fn((n_samples, n_features), |_| rng.gen_range(0..n_categories))
+}
+
+fn k_modes_bench(c: &mut Criterion) {
+    let mut rng = Xoshiro256Plus::seed_from_u64(42);
+    let sample_sizes = [1_000, 10_000, 20_000];
+    let feature_dims = [5, 10];
+    let n_clusters = 4;
+    let n_categories = 5;
+
+    let mut benchmark = c.benchmark_group("k_modes");
+    config::set_default_benchmark_configs(&mut benchmark);
+    benchmark.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    for &n_features in &feature_dims {
+        for &n_samples in &sample_sizes {
+            let raw_data =
+                generate_synthetic_categorical(n_samples, n_features, n_categories, &mut rng);
+            let dataset = DatasetBase::from(raw_data);
+
+            // Benchmark Cao initialization
+            benchmark.bench_function(
+                BenchmarkId::new("Cao", format!("{n_samples}x{n_features}")),
+                |bencher| {
+                    bencher.iter(|| {
+                        black_box(
+                            KModes::params(black_box(n_clusters))
+                                .max_n_iterations(black_box(20))
+                                .init_method(KModesInit::Cao)
+                                .fit(&dataset)
+                                .unwrap(),
+                        )
+                    });
+                },
+            );
+
+            // Benchmark Huang initialization
+            benchmark.bench_function(
+                BenchmarkId::new("Huang", format!("{n_samples}x{n_features}")),
+                |bencher| {
+                    bencher.iter(|| {
+                        black_box(
+                            KModes::params(black_box(n_clusters))
+                                .max_n_iterations(black_box(20))
+                                .init_method(KModesInit::Huang)
+                                .n_runs(1)
+                                .fit(&dataset)
+                                .unwrap(),
+                        )
+                    });
+                },
+            );
+
+            // Benchmark Random initialization
+            benchmark.bench_function(
+                BenchmarkId::new("Random", format!("{n_samples}x{n_features}")),
+                |bencher| {
+                    bencher.iter(|| {
+                        black_box(
+                            KModes::params(black_box(n_clusters))
+                                .max_n_iterations(black_box(20))
+                                .init_method(KModesInit::Random)
+                                .n_runs(1)
+                                .fit(&dataset)
+                                .unwrap(),
+                        )
+                    });
+                },
+            );
+        }
+    }
+
+    benchmark.finish();
+}
+
+#[cfg(not(target_os = "windows"))]
+criterion_group! {
+    name = benches;
+    config = config::get_default_profiling_configs();
+    targets = k_modes_bench
+}
+#[cfg(target_os = "windows")]
+criterion_group!(benches, k_modes_bench);
+
+criterion_main!(benches);
diff --git a/algorithms/linfa-clustering/src/k_modes/algorithm.rs b/algorithms/linfa-clustering/src/k_modes/algorithm.rs
new file mode 100644
index 000000000..504e96d67
--- /dev/null
+++ b/algorithms/linfa-clustering/src/k_modes/algorithm.rs
@@ -0,0 +1,964 @@
+use std::collections::HashMap;
+
+use crate::k_modes::init::KModesInit;
+use crate::k_modes::{KModesError, KModesParams, KModesValidParams};
+use linfa::traits::{Fit, PredictInplace};
+use linfa::DatasetBase;
+use ndarray::{Array1, Array2, ArrayBase, ArrayView1, ArrayView2, Data, Ix1, Ix2};
+use ndarray_rand::rand::Rng;
+use rand_xoshiro::Xoshiro256Plus;
+#[cfg(feature = "serde")]
+use serde_crate::{Deserialize, Serialize};
+
+#[cfg_attr(
+    feature = "serde",
+    derive(Serialize, Deserialize),
+    serde(crate = "serde_crate")
+)]
+#[derive(Clone, Debug, PartialEq, Eq)]
+/// K-Modes clustering model for categorical data (Huang \[1998\]).
+///
+/// Partitions $N$ categorical observations into $K$ clusters by finding a mode vector $Q_l$
+/// for each cluster $l \in \{0, \dots, K-1\}$ that minimizes the sum of matching dissimilarities:
+///
+/// $$P(W, Q) = \sum_{l=0}^{K-1} \sum_{i=1}^{N} w_{il} \sum_{j=1}^{M} \delta(x_{ij}, q_{lj})$$
+///
+/// where:
+/// - $w_{il} \in \{0, 1\}$ indicates whether observation $X_i$ belongs to cluster $l$,
+/// - $\delta(x_{ij}, q_{lj}) = 0$ if $x_{ij} = q_{lj}$, and $1$ otherwise (Simple Matching Dissimilarity),
+/// - $q_{lj}$ is the empirical mode (most frequent category) of attribute $j$ in cluster $l$.
+///
+/// ## Tutorial
+///
+/// ```
+/// use linfa::DatasetBase;
+/// use linfa::traits::{Fit, Predict};
+/// use linfa_clustering::{KModes, KModesInit};
+/// use ndarray::array;
+///
+/// // Create a categorical dataset
+/// let entries = array![
+///     ["Sunny", "Hot", "High"],
+///     ["Sunny", "Hot", "High"],
+///     ["Rainy", "Mild", "High"],
+///     ["Rainy", "Cool", "Normal"],
+///     ["Rainy", "Cool", "Normal"],
+///     ["Sunny", "Hot", "High"]
+/// ];
+/// let dataset = DatasetBase::from(entries);
+///
+/// // Configure and fit K-Modes with 2 clusters
+/// let model = KModes::params(2)
+///     .max_n_iterations(100)
+///     .init_method(KModesInit::Cao)
+///     .fit(&dataset)
+///     .expect("K-Modes fitting failed");
+///
+/// // Predict cluster assignments
+/// let predictions = model.predict(&dataset);
+/// assert_eq!(predictions.len(), 6);
+///
+/// // Query fitted cluster modes and total cost
+/// assert_eq!(model.modes().dim(), (2, 3));
+/// assert!(model.cost() <= 6);
+/// ```
+pub struct KModes<T: EquivalenceTarget> {
+    modes: Array2<T>,
+    cost: usize,
+}
+
+/// Trait bound for categorical elements supported by K-Modes.
+pub trait EquivalenceTarget: PartialEq + Eq + Clone + std::hash::Hash {}
+impl<T: PartialEq + Eq + Clone + std::hash::Hash> EquivalenceTarget for T {}
+
+impl<T: EquivalenceTarget> KModes<T> {
+    /// Constructs a parameter builder with `n_clusters` clusters and default RNG.
+    pub fn params(n_clusters: usize) -> KModesParams<T, Xoshiro256Plus> {
+        KModesParams::new(n_clusters)
+    }
+
+    /// Constructs a parameter builder with `n_clusters` clusters and a custom RNG.
+    pub fn params_with_rng<R: Rng>(n_clusters: usize, rng: R) -> KModesParams<T, R> {
+        KModesParams::new_with_rng(n_clusters, rng)
+    }
+
+    /// Cluster modes matrix with shape `(n_clusters, n_features)`.
+    pub fn modes(&self) -> &Array2<T> {
+        &self.modes
+    }
+
+    /// Total dissimilarity cost of the fitted clustering.
+    pub fn cost(&self) -> usize {
+        self.cost
+    }
+}
+
+impl<T: EquivalenceTarget, R: Rng + Clone, D: Data<Elem = T>, L>
+    Fit<ArrayBase<D, Ix2>, L, KModesError> for KModesValidParams<T, R>
+{
+    type Object = KModes<T>;
+
+    /// Fits the K-Modes model on a categorical dataset.
+    ///
+    /// Executes `n_runs` independent clusterings (or 1 run for deterministic `Cao` initialization)
+    /// and returns the model with minimal total dissimilarity cost.
+    fn fit(
+        &self,
+        dataset: &DatasetBase<ArrayBase<D, Ix2>, L>,
+    ) -> Result<Self::Object, KModesError> {
+        let observations = dataset.records().view();
+        let (n_points, _) = observations.dim();
+
+        if n_points == 0 {
+            return Err(
+                linfa::error::Error::Parameters("Dataset cannot be empty".to_string()).into(),
+            );
+        }
+
+        if self.n_clusters() > n_points {
+            return Err(linfa::error::Error::Parameters(format!(
+                "Cannot have more clusters ({}) than data points ({})",
+                self.n_clusters(),
+                n_points
+            ))
+            .into());
+        }
+
+        let mut rng = self.rng().clone();
+        let mut best_cost = usize::MAX;
+        let mut best_centroids = None;
+
+        let n_runs = match self.init_method() {
+            KModesInit::Cao if self.n_runs() > 1 => {
+                if self.verbose() {
+                    println!("Cao initialization is deterministic. Running 1 initialization.");
+                }
+                1
+            }
+            KModesInit::Precomputed(_) if self.n_runs() > 1 => 1,
+            _ => self.n_runs(),
+        };
+
+        for run_idx in 0..n_runs {
+            if self.verbose() && n_runs > 1 {
+                println!("Starting K-Modes run {}/{}", run_idx + 1, n_runs);
+            }
+
+            let (centroids, cost) = k_modes_single(
+                observations,
+                self.n_clusters(),
+                self.max_n_iterations(),
+                self.init_method(),
+                &mut rng,
+                self.verbose(),
+            );
+
+            if cost < best_cost {
+                best_cost = cost;
+                best_centroids = Some(centroids);
+            }
+        }
+
+        let modes = best_centroids.expect("internal error: K-Modes failed to fit centroids");
+
+        Ok(KModes {
+            modes,
+            cost: best_cost,
+        })
+    }
+}
+
+/// Simple matching dissimilarity: counts mismatched features between two categorical vectors.
+///
+/// $$d(X, Y) = \sum_{j=1}^{M} [x_j \neq y_j]$$
+pub(crate) fn matching_dissim<T: PartialEq>(a: ArrayView1<T>, b: ArrayView1<T>) -> usize {
+    a.iter().zip(b.iter()).filter(|(x, y)| x != y).count()
+}
+
+/// Finds the nearest centroid index and minimum distance for a given observation point.
+pub(crate) fn closest_centroid<T: PartialEq>(
+    centroids: ArrayView2<T>,
+    observation: ArrayView1<T>,
+) -> (usize, usize) {
+    let mut min_dist = usize::MAX;
+    let mut closest_idx = 0;
+
+    for (idx, centroid) in centroids.rows().into_iter().enumerate() {
+        let dist = matching_dissim(centroid, observation);
+        if dist < min_dist {
+            min_dist = dist;
+            closest_idx = idx;
+        }
+    }
+
+    (closest_idx, min_dist)
+}
+
+/// Evaluates total clustering loss: sum of matching distances of all observations to their assigned modes.
+fn calculate_cost<T: PartialEq>(
+    observations: ArrayView2<T>,
+    centroids: ArrayView2<T>,
+    memberships: &[usize],
+) -> usize {
+    observations
+        .rows()
+        .into_iter()
+        .zip(memberships.iter())
+        .map(|(obs, &cluster_id)| matching_dissim(obs, centroids.row(cluster_id)))
+        .sum()
+}
+
+/// Single run of Huang's incremental online K-Modes algorithm.
+///
+/// Steps:
+/// 1. Initialize $K$ modes using the configured initialization strategy (`Cao`, `Huang`, `Random`, `Precomputed`).
+/// 2. Assign each observation to its closest initial mode and build frequency tables:
+///    `cl_attr_freq[cluster][attribute][value] -> count`
+/// 3. Recompute initial cluster modes by taking the argmax frequency per attribute.
+/// 4. Iteratively cycle through all points:
+///    - Determine the closest cluster mode for point $X_i$.
+///    - If closest cluster $\neq$ current cluster, immediately move $X_i$, updating the category
+///      frequency counts and modes for both the source and target clusters in $O(M)$ time.
+///    - If the source cluster becomes empty, reinitialize it with a point drawn from the largest cluster.
+/// 5. Repeat until no points change cluster membership (`moves == 0`) or `max_n_iterations` is reached.
+fn k_modes_single<T: EquivalenceTarget, R: Rng>(
+    observations: ArrayView2<T>,
+    n_clusters: usize,
+    max_n_iterations: u64,
+    init: &KModesInit<T>,
+    rng: &mut R,
+    verbose: bool,
+) -> (Array2<T>, usize) {
+    let (n_points, n_attrs) = observations.dim();
+
+    // centroid generation
+    let mut centroids = init.run(n_clusters, observations, rng);
+    let mut memberships = vec![0usize; n_points];
+    let mut cluster_counts = vec![0usize; n_clusters];
+
+    // Track category frequencies per cluster and feature: cl_attr_freq[k][j][v] -> count
+    let mut cl_attr_freq: Vec<Vec<HashMap<T, usize>>> = (0..n_clusters)
+        .map(|_| (0..n_attrs).map(|_| HashMap::new()).collect())
+        .collect();
+
+    for (ipoint, curpoint) in observations.rows().into_iter().enumerate() {
+        let (clust, _) = closest_centroid(centroids.view(), curpoint);
+        memberships[ipoint] = clust;
+        cluster_counts[clust] += 1;
+
+        for (iattr, curattr) in curpoint.iter().enumerate() {
+            *cl_attr_freq[clust][iattr]
+                .entry(curattr.clone())
+                .or_insert(0) += 1;
+        }
+    }
+
+    // Step 3: Compute initial empirical modes from frequency distributions
+    for ik in 0..n_clusters {
+        if cluster_counts[ik] == 0 {
+            let random_idx = rng.gen_range(0..n_points);
+            let sample_row = observations.row(random_idx);
+            for iattr in 0..n_attrs {
+                centroids[[ik, iattr]] = sample_row[iattr].clone();
+            }
+        } else {
+            for iattr in 0..n_attrs {
+                let cur_val = &centroids[[ik, iattr]];
+                let cur_count = cl_attr_freq[ik][iattr].get(cur_val).copied().unwrap_or(0);
+                if let Some((best_val, &max_count)) = cl_attr_freq[ik][iattr]
+                    .iter()
+                    .filter(|(_, &count)| count > 0)
+                    .max_by_key(|(_, &count)| count)
+                {
+                    if max_count > cur_count {
+                        centroids[[ik, iattr]] = best_val.clone();
+                    }
+                }
+            }
+        }
+    }
+
+    let mut cost = calculate_cost(observations, centroids.view(), &memberships);
+
+    // Online incremental reassignment loop
+    let mut iter = 0;
+    while iter < max_n_iterations {
+        iter += 1;
+        let mut moves = 0;
+
+        for ipoint in 0..n_points {
+            let curpoint = observations.row(ipoint);
+            let (to_clust, _) = closest_centroid(centroids.view(), curpoint);
+            let from_clust = memberships[ipoint];
+
+            if to_clust == from_clust {
+                continue;
+            }
+
+            // Move point to closer cluster
+            moves += 1;
+            memberships[ipoint] = to_clust;
+            cluster_counts[from_clust] -= 1;
+            cluster_counts[to_clust] += 1;
+
+            // Dynamically update frequency maps and centroid modes
+            move_point_cat(
+                curpoint,
+                to_clust,
+                from_clust,
+                &mut cl_attr_freq,
+                &mut centroids,
+            );
+
+            // Reassign a point from the largest cluster if a cluster became empty
+            if cluster_counts[from_clust] == 0 {
+                let largest_clust = cluster_counts
+                    .iter()
+                    .enumerate()
+                    .max_by_key(|(_, &count)| count)
+                    .map(|(idx, _)| idx)
+                    .unwrap_or(0);
+
+                let candidates: Vec<usize> = memberships
+                    .iter()
+                    .enumerate()
+                    .filter_map(|(idx, &c)| if c == largest_clust { Some(idx) } else { None })
+                    .collect();
+
+                if !candidates.is_empty() {
+                    let rindx = candidates[rng.gen_range(0..candidates.len())];
+                    let rpoint = observations.row(rindx);
+
+                    memberships[rindx] = from_clust;
+                    cluster_counts[largest_clust] -= 1;
+                    cluster_counts[from_clust] += 1;
+
+                    move_point_cat(
+                        rpoint,
+                        from_clust,
+                        largest_clust,
+                        &mut cl_attr_freq,
+                        &mut centroids,
+                    );
+                }
+            }
+        }
+
+        let ncost = calculate_cost(observations, centroids.view(), &memberships);
+        cost = ncost;
+
+        // Convergence check: loop terminates when no cluster reassignments occur during the epoch
+        if moves == 0 {
+            if verbose {
+                println!(
+                    "K-Modes converged at iteration {} (moves: {}, cost: {})",
+                    iter, moves, cost
+                );
+            }
+            break;
+        } else if verbose && (iter % 10 == 0 || iter == max_n_iterations) {
+            println!(
+                "Iteration {}/{}: moves = {}, cost = {}",
+                iter, max_n_iterations, moves, cost
+            );
+        }
+    }
+
+    (centroids, cost)
+}
+
+/// Dynamically updates attribute frequencies and modes when transferring a single observation between clusters.
+///
+/// - For `to_clust`: increments category frequency; if new frequency exceeds the current mode's count,
+///   the centroid coordinate is updated to this category in $O(1)$ time.
+/// - For `from_clust`: decrements category frequency; if the removed category was the current mode,
+///   the new mode is recalculated by finding the key with maximum count.
+fn move_point_cat<T: EquivalenceTarget>(
+    point: ArrayView1<T>,
+    to_clust: usize,
+    from_clust: usize,
+    cl_attr_freq: &mut [Vec<HashMap<T, usize>>],
+    centroids: &mut Array2<T>,
+) {
+    for (iattr, curattr) in point.iter().enumerate() {
+        // Increment frequency in destination cluster
+        let to_map = &mut cl_attr_freq[to_clust][iattr];
+        let to_count = to_map.entry(curattr.clone()).or_insert(0);
+        *to_count += 1;
+        let new_val_freq = *to_count;
+
+        let current_centroid_val = &centroids[[to_clust, iattr]];
+        let current_centroid_freq = to_map.get(current_centroid_val).copied().unwrap_or(0);
+        if new_val_freq > current_centroid_freq {
+            centroids[[to_clust, iattr]] = curattr.clone();
+        }
+
+        // Decrement frequency in source cluster
+        let from_map = &mut cl_attr_freq[from_clust][iattr];
+        if let Some(from_count) = from_map.get_mut(curattr) {
+            if *from_count > 0 {
+                *from_count -= 1;
+            }
+        }
+
+        // If the departed value was the mode, check if another value is now strictly more frequent
+        let old_centroid_val = &centroids[[from_clust, iattr]];
+        if old_centroid_val == curattr {
+            let remaining_count = from_map.get(curattr).copied().unwrap_or(0);
+            if let Some((best_val, &max_count)) = from_map
+                .iter()
+                .filter(|(_, &count)| count > 0)
+                .max_by_key(|(_, &count)| count)
+            {
+                if max_count > remaining_count {
+                    centroids[[from_clust, iattr]] = best_val.clone();
+                }
+            }
+        }
+    }
+}
+
+impl<T: EquivalenceTarget, D: Data<Elem = T>> PredictInplace<ArrayBase<D, Ix2>, Array1<usize>>
+    for KModes<T>
+{
+    /// Predicts closest cluster indices for a 2D batch of observations.
+    fn predict_inplace(&self, observations: &ArrayBase<D, Ix2>, memberships: &mut Array1<usize>) {
+        assert_eq!(observations.nrows(), memberships.len());
+
+        for (i, obs) in observations.rows().into_iter().enumerate() {
+            memberships[i] = closest_centroid(self.modes.view(), obs).0;
+        }
+    }
+
+    fn default_target(&self, x: &ArrayBase<D, Ix2>) -> Array1<usize> {
+        Array1::zeros(x.nrows())
+    }
+}
+
+impl<T: EquivalenceTarget, D: Data<Elem = T>> PredictInplace<ArrayBase<D, Ix1>, usize>
+    for KModes<T>
+{
+    /// Predicts closest cluster index for a single 1D observation.
+    fn predict_inplace(&self, observation: &ArrayBase<D, Ix1>, membership: &mut usize) {
+        *membership = closest_centroid(self.modes.view(), observation.view()).0;
+    }
+
+    fn default_target(&self, _x: &ArrayBase<D, Ix1>) -> usize {
+        0
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use linfa::traits::Predict;
+    use ndarray::array;
+    use ndarray_rand::rand::SeedableRng;
+
+    #[test]
+    fn autotraits() {
+        fn has_autotraits<T: Send + Sync + Sized + Unpin>() {}
+        has_autotraits::<KModes<usize>>();
+        has_autotraits::<KModes<&str>>();
+        has_autotraits::<KModes<String>>();
+        has_autotraits::<KModes<u8>>();
+        has_autotraits::<KModes<char>>();
+        has_autotraits::<KModes<bool>>();
+    }
+
+    #[test]
+    fn test_kmodes_simple_fitting() {
+        let data = array![
+            ["A", "X"],
+            ["A", "X"],
+            ["B", "Y"],
+            ["B", "Y"],
+            ["A", "X"],
+            ["B", "Y"]
+        ];
+
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2)
+            .max_n_iterations(50)
+            .verbose(false)
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes().dim(), (2, 2));
+
+        let predictions = model.predict(&dataset);
+        assert_eq!(predictions.len(), 6);
+
+        assert_eq!(predictions[0], predictions[1]);
+        assert_eq!(predictions[0], predictions[4]);
+
+        assert_eq!(predictions[2], predictions[3]);
+        assert_eq!(predictions[2], predictions[5]);
+        assert_ne!(predictions[0], predictions[2]);
+        assert_eq!(model.cost(), 0);
+    }
+
+    #[test]
+    fn test_kmodes_single_observation_predict() {
+        let data = array![["A", "X"], ["B", "Y"]];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        let single_sample = array!["A", "X"];
+        let pred: usize = model.predict(&single_sample);
+        assert!(pred < 2);
+    }
+
+    #[test]
+    fn test_kmodes_empty_dataset_error() {
+        let empty_data: Array2<usize> = Array2::zeros((0, 3));
+        let dataset = DatasetBase::from(empty_data);
+        let result = KModes::params(2).fit(&dataset);
+        assert!(result.is_err());
+    }
+
+    #[test]
+    fn test_kmodes_n_clusters_exceeds_samples_error() {
+        let data = array![["A", "X"], ["B", "Y"]];
+        let dataset = DatasetBase::from(data);
+        let result = KModes::params(5).fit(&dataset);
+        assert!(result.is_err());
+    }
+
+    #[test]
+    fn test_kmodes_single_cluster() {
+        let data = array![["A", "X"], ["A", "X"], ["A", "Y"], ["B", "X"]];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(1).fit(&dataset).unwrap();
+
+        assert_eq!(model.modes().dim(), (1, 2));
+        assert_eq!(model.modes()[[0, 0]], "A");
+        assert_eq!(model.modes()[[0, 1]], "X");
+
+        let predictions = model.predict(&dataset);
+        assert_eq!(predictions.to_vec(), vec![0, 0, 0, 0]);
+    }
+
+    #[test]
+    fn test_kmodes_n_clusters_equals_n_samples() {
+        let data = array![["A", "X"], ["B", "Y"], ["C", "Z"]];
+        let dataset = DatasetBase::from(data.clone());
+        let model = KModes::params(3)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes().dim(), (3, 2));
+        assert_eq!(model.cost(), 0);
+
+        let predictions = model.predict(&dataset);
+        let mut unique_labels: Vec<usize> = predictions.to_vec();
+        unique_labels.sort_unstable();
+        unique_labels.dedup();
+        assert_eq!(unique_labels.len(), 3);
+    }
+
+    #[test]
+    fn test_kmodes_identical_observations() {
+        let data = array![
+            ["A", "X", "1"],
+            ["A", "X", "1"],
+            ["A", "X", "1"],
+            ["A", "X", "1"]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        assert_eq!(model.cost(), 0);
+        let predictions = model.predict(&dataset);
+        assert_eq!(predictions.len(), 4);
+    }
+
+    #[test]
+    fn test_kmodes_high_dimensional_dataset() {
+        let row_a = vec!["A"; 20];
+        let row_b = vec!["B"; 20];
+        let mut raw_data = Vec::new();
+        for _ in 0..10 {
+            raw_data.extend(row_a.clone());
+        }
+        for _ in 0..10 {
+            raw_data.extend(row_b.clone());
+        }
+        let data = Array2::from_shape_vec((20, 20), raw_data).unwrap();
+        let dataset = DatasetBase::from(data);
+
+        let model = KModes::params(2)
+            .max_n_iterations(20)
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes().dim(), (2, 20));
+        assert_eq!(model.cost(), 0);
+
+        let preds = model.predict(&dataset);
+        assert_eq!(preds.slice(ndarray::s![0..10]).to_vec(), vec![preds[0]; 10]);
+        assert_eq!(
+            preds.slice(ndarray::s![10..20]).to_vec(),
+            vec![preds[10]; 10]
+        );
+        assert_ne!(preds[0], preds[10]);
+    }
+
+    #[test]
+    fn test_kmodes_max_iterations_bound() {
+        let data = array![
+            ["A", "1"],
+            ["A", "2"],
+            ["B", "1"],
+            ["B", "2"],
+            ["C", "1"],
+            ["C", "2"]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(3).max_n_iterations(1).fit(&dataset).unwrap();
+
+        assert_eq!(model.modes().dim(), (3, 2));
+    }
+
+    #[test]
+    fn test_kmodes_n_runs_selects_minimal_cost() {
+        let data = array![
+            ["A", "X", "1"],
+            ["A", "X", "2"],
+            ["B", "Y", "1"],
+            ["B", "Y", "2"],
+            ["C", "Z", "1"],
+            ["C", "Z", "2"]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model_multi = KModes::params(3)
+            .init_method(KModesInit::Random)
+            .n_runs(10)
+            .fit(&dataset)
+            .unwrap();
+
+        let model_single = KModes::params(3)
+            .init_method(KModesInit::Random)
+            .n_runs(1)
+            .fit(&dataset)
+            .unwrap();
+
+        assert!(model_multi.cost() <= model_single.cost() || model_single.cost() == 0);
+    }
+
+    #[test]
+    fn test_kmodes_with_custom_rng() {
+        let data = array![
+            ["A", "X"],
+            ["A", "X"],
+            ["B", "Y"],
+            ["B", "Y"],
+            ["C", "Z"],
+            ["C", "Z"]
+        ];
+        let dataset = DatasetBase::from(data);
+
+        let rng1 = Xoshiro256Plus::seed_from_u64(42);
+        let rng2 = Xoshiro256Plus::seed_from_u64(42);
+
+        let model1 = KModes::params_with_rng(3, rng1)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+
+        let model2 = KModes::params_with_rng(3, rng2)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model1.modes(), model2.modes());
+        assert_eq!(model1.cost(), model2.cost());
+    }
+
+    #[test]
+    fn test_kmodes_predict_unseen_categories() {
+        let data = array![["A", "X"], ["A", "X"], ["B", "Y"], ["B", "Y"]];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        // Sample with category "Z" never seen in training
+        let unseen = array![["A", "Z"]];
+        let pred: usize = model.predict(&unseen.row(0));
+        let expected_cluster = model.predict(&array!["A", "X"]);
+        assert_eq!(pred, expected_cluster);
+    }
+
+    #[test]
+    fn test_kmodes_types_numeric() {
+        let data = array![
+            [1usize, 10usize],
+            [1usize, 10usize],
+            [2usize, 20usize],
+            [2usize, 20usize]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2).fit(&dataset).unwrap();
+        assert_eq!(model.modes().dim(), (2, 2));
+        assert_eq!(model.cost(), 0);
+
+        let u8_data = array![[1u8, 2u8], [1u8, 2u8], [3u8, 4u8], [3u8, 4u8]];
+        let u8_dataset = DatasetBase::from(u8_data);
+        let u8_model = KModes::params(2).fit(&u8_dataset).unwrap();
+        assert_eq!(u8_model.modes().dim(), (2, 2));
+        assert_eq!(u8_model.cost(), 0);
+
+        let i32_data = array![
+            [-10i32, 100i32],
+            [-10i32, 100i32],
+            [50i32, -200i32],
+            [50i32, -200i32]
+        ];
+        let i32_dataset = DatasetBase::from(i32_data);
+        let i32_model = KModes::params(2).fit(&i32_dataset).unwrap();
+        assert_eq!(i32_model.modes().dim(), (2, 2));
+        assert_eq!(i32_model.cost(), 0);
+    }
+
+    #[test]
+    fn test_kmodes_types_char_and_bool() {
+        let char_data = array![['a', 'x'], ['a', 'x'], ['b', 'y'], ['b', 'y']];
+        let char_dataset = DatasetBase::from(char_data);
+        let char_model = KModes::params(2).fit(&char_dataset).unwrap();
+        assert_eq!(char_model.modes().dim(), (2, 2));
+        assert_eq!(char_model.cost(), 0);
+
+        let bool_data = array![[true, false], [true, false], [false, true], [false, true]];
+        let bool_dataset = DatasetBase::from(bool_data);
+        let bool_model = KModes::params(2).fit(&bool_dataset).unwrap();
+        assert_eq!(bool_model.modes().dim(), (2, 2));
+        assert_eq!(bool_model.cost(), 0);
+    }
+
+    #[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)]
+    enum Weather {
+        Sunny,
+        Rainy,
+        Overcast,
+    }
+
+    #[test]
+    fn test_kmodes_types_custom_enum() {
+        let data = array![
+            [Weather::Sunny, Weather::Sunny],
+            [Weather::Sunny, Weather::Sunny],
+            [Weather::Rainy, Weather::Rainy],
+            [Weather::Rainy, Weather::Rainy],
+            [Weather::Overcast, Weather::Overcast],
+            [Weather::Overcast, Weather::Overcast]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(3)
+            .init_method(KModesInit::Cao)
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes().dim(), (3, 2));
+        assert_eq!(model.cost(), 0);
+
+        let preds = model.predict(&dataset);
+        assert_eq!(preds[0], preds[1]);
+        assert_eq!(preds[2], preds[3]);
+        assert_eq!(preds[4], preds[5]);
+    }
+
+    #[test]
+    fn test_kmodes_predict_dataset_and_records_consistency() {
+        let data = array![
+            ["Cat", "Small"],
+            ["Dog", "Large"],
+            ["Cat", "Small"],
+            ["Dog", "Large"]
+        ];
+        let dataset = DatasetBase::from(data.clone());
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        let dataset_preds = model.predict(&dataset);
+        let records_preds = model.predict(&data);
+
+        assert_eq!(dataset_preds, records_preds);
+        for (i, row) in data.rows().into_iter().enumerate() {
+            let single_pred: usize = model.predict(&row);
+            assert_eq!(dataset_preds[i], single_pred);
+        }
+    }
+
+    #[test]
+    fn test_kmodes_all_init_strategies() {
+        let data = array![
+            ["A", "1"],
+            ["A", "1"],
+            ["B", "2"],
+            ["B", "2"],
+            ["C", "3"],
+            ["C", "3"]
+        ];
+        let dataset = DatasetBase::from(data);
+
+        let m_cao = KModes::params(3)
+            .init_method(KModesInit::Cao)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_cao.modes().dim(), (3, 2));
+
+        let m_huang = KModes::params(3)
+            .init_method(KModesInit::Huang)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_huang.modes().dim(), (3, 2));
+
+        let m_random = KModes::params(3)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_random.modes().dim(), (3, 2));
+
+        let precomputed = array![["A", "1"], ["B", "2"], ["C", "3"]];
+        let m_pre = KModes::params(3)
+            .init_method(KModesInit::Precomputed(precomputed))
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_pre.modes().dim(), (3, 2));
+        assert_eq!(m_pre.cost(), 0);
+    }
+
+    #[cfg(feature = "serde")]
+    #[test]
+    fn test_kmodes_serde() {
+        let data = array![
+            ["A".to_string(), "X".to_string()],
+            ["A".to_string(), "X".to_string()],
+            ["B".to_string(), "Y".to_string()],
+            ["B".to_string(), "Y".to_string()]
+        ];
+        let dataset = DatasetBase::from(data);
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        let serialized = serde_json::to_string(&model).unwrap();
+        let deserialized: KModes<String> = serde_json::from_str(&serialized).unwrap();
+
+        assert_eq!(model, deserialized);
+        assert_eq!(model.modes(), deserialized.modes());
+        assert_eq!(model.cost(), deserialized.cost());
+    }
+
+    #[test]
+    fn test_kmodes_verbose_logging_coverage() {
+        let data = array![
+            ["A", "X"],
+            ["A", "X"],
+            ["B", "Y"],
+            ["B", "Y"],
+            ["C", "Z"],
+            ["C", "Z"]
+        ];
+        let dataset = DatasetBase::from(data);
+
+        // Triggers Cao deterministic warning logging
+        let m_cao = KModes::params(2)
+            .n_runs(3)
+            .verbose(true)
+            .init_method(KModesInit::Cao)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_cao.modes().dim(), (2, 2));
+
+        // Triggers multi-run and convergence verbose logging
+        let m_rand = KModes::params(2)
+            .n_runs(2)
+            .max_n_iterations(20)
+            .verbose(true)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_rand.modes().dim(), (2, 2));
+
+        // Triggers max_n_iterations logging branch
+        let m_max_iter = KModes::params(2)
+            .n_runs(1)
+            .max_n_iterations(1)
+            .verbose(true)
+            .init_method(KModesInit::Random)
+            .fit(&dataset)
+            .unwrap();
+        assert_eq!(m_max_iter.modes().dim(), (2, 2));
+    }
+
+    #[test]
+    fn test_kmodes_empty_cluster_recovery() {
+        let data = array![
+            ["A", "1"],
+            ["A", "1"],
+            ["A", "1"],
+            ["A", "1"],
+            ["B", "2"],
+            ["B", "2"],
+            ["B", "2"],
+            ["B", "2"]
+        ];
+        let dataset = DatasetBase::from(data);
+
+        // Precompute centroids with duplicate points so cluster 1 becomes empty during assignment
+        let init_modes = array![["A", "1"], ["A", "1"], ["B", "2"]];
+        let model = KModes::params(3)
+            .init_method(KModesInit::Precomputed(init_modes))
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes().dim(), (3, 2));
+    }
+
+    #[test]
+    fn test_kmodes_default_target() {
+        use linfa::traits::PredictInplace;
+        let data = array![["A", "X"], ["B", "Y"]];
+        let dataset = DatasetBase::from(data.clone());
+        let model = KModes::params(2).fit(&dataset).unwrap();
+
+        let t2d = model.default_target(&data);
+        assert_eq!(t2d.len(), 2);
+
+        let sample = array!["A", "X"];
+        let t1d = model.default_target(&sample);
+        assert_eq!(t1d, 0);
+    }
+
+    #[test]
+    fn test_kmodes_initial_empirical_mode_update() {
+        // Points assigned to cluster 0 will have "Y" as mode, differing from initial seed "X"
+        let data = array![
+            ["Y", "1"],
+            ["Y", "1"],
+            ["Y", "1"],
+            ["X", "1"],
+            ["Z", "9"],
+            ["Z", "9"]
+        ];
+        let dataset = DatasetBase::from(data);
+
+        // Precomputed initial centroids: cluster 0 starts with ["X", "1"], cluster 1 starts with ["Z", "9"]
+        let init_modes = array![["X", "1"], ["Z", "9"]];
+        let model = KModes::params(2)
+            .init_method(KModesInit::Precomputed(init_modes))
+            .fit(&dataset)
+            .unwrap();
+
+        assert_eq!(model.modes()[[0, 0]], "Y");
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_kmodes_predict_inplace_mismatched_length() {
+        use linfa::traits::PredictInplace;
+        let data = array![["A", "1"], ["B", "2"]];
+        let dataset = DatasetBase::from(data.clone());
+        let model = KModes::params(2).fit(&dataset).unwrap();
+        let mut wrong_memberships = ndarray::Array1::zeros(10);
+        model.predict_inplace(&data, &mut wrong_memberships);
+    }
+}
diff --git a/algorithms/linfa-clustering/src/k_modes/errors.rs b/algorithms/linfa-clustering/src/k_modes/errors.rs
new file mode 100644
index 000000000..dda92c7d0
--- /dev/null
+++ b/algorithms/linfa-clustering/src/k_modes/errors.rs
@@ -0,0 +1,19 @@
+use thiserror::Error;
+
+#[derive(Error, Debug, PartialEq)]
+pub enum KModesParamsError {
+    #[error("n_clusters cannot be 0")]
+    NClusters,
+    #[error("max_n_iterations cannot be 0")]
+    MaxIterations,
+    #[error("n_runs cannot be 0")]
+    NRuns,
+}
+
+#[derive(Error, Debug)]
+pub enum KModesError {
+    #[error("Invalid hyperparameter: {0}")]
+    InvalidParams(#[from] KModesParamsError),
+    #[error(transparent)]
+    LinfaError(#[from] linfa::error::Error),
+}
diff --git a/algorithms/linfa-clustering/src/k_modes/hyperparams.rs b/algorithms/linfa-clustering/src/k_modes/hyperparams.rs
new file mode 100644
index 000000000..cf4a5cafb
--- /dev/null
+++ b/algorithms/linfa-clustering/src/k_modes/hyperparams.rs
@@ -0,0 +1,273 @@
+use crate::k_modes::init::KModesInit;
+use crate::KModesParamsError;
+use linfa::ParamGuard;
+use ndarray_rand::rand::{Rng, SeedableRng};
+use rand_xoshiro::Xoshiro256Plus;
+#[cfg(feature = "serde")]
+use serde_crate::{Deserialize, Serialize};
+
+#[cfg_attr(
+    feature = "serde",
+    derive(Serialize, Deserialize),
+    serde(crate = "serde_crate")
+)]
+#[derive(Clone, Debug, PartialEq)]
+/// Hyperparameters that can be specified for
+/// the [K-Modes algorithm](crate::KModes).
+pub struct KModesValidParams<T, R: Rng> {
+    /// Maximum number of iterations for a single run
+    /// When max_n_iterations is exceeded we terminate training
+    pub(crate) max_n_iterations: u64,
+    /// Number of clusters to form
+    pub(crate) n_clusters: usize,
+    /// Number of times the algorithm will run with different seeds
+    pub(crate) n_runs: usize,
+    /// Centroid initialization methods
+    pub(crate) init: KModesInit<T>,
+    /// Enable verbose progress logging to stdout.
+    pub(crate) verbose: bool,
+    /// Random number generator.
+    pub(crate) rng: R,
+}
+
+impl<T, R: Rng> KModesValidParams<T, R> {
+    pub fn n_clusters(&self) -> usize {
+        self.n_clusters
+    }
+
+    pub fn max_n_iterations(&self) -> u64 {
+        self.max_n_iterations
+    }
+
+    pub fn n_runs(&self) -> usize {
+        self.n_runs
+    }
+
+    pub fn init_method(&self) -> &KModesInit<T> {
+        &self.init
+    }
+
+    pub fn verbose(&self) -> bool {
+        self.verbose
+    }
+
+    pub fn rng(&self) -> &R {
+        &self.rng
+    }
+}
+
+#[cfg_attr(
+    feature = "serde",
+    derive(Serialize, Deserialize),
+    serde(crate = "serde_crate")
+)]
+#[derive(Clone, Debug, PartialEq)]
+/// Helper builder to configure hyperparameters for the [K-Modes algorithm](crate::KModes).
+pub struct KModesParams<T, R: Rng>(KModesValidParams<T, R>);
+
+impl<T> KModesParams<T, Xoshiro256Plus> {
+    /// Create a new K-Modes parameter builder with default RNG (seed 42).
+    /// Implemented with `new_with_rng` with default parameters:
+    ///  * max_n_iterations: 100
+    ///  * n_runs: 10
+    ///  * init: KModesInit::Cao
+    ///  * verbose: false
+    pub fn new(n_clusters: usize) -> Self {
+        Self::new_with_rng(n_clusters, Xoshiro256Plus::seed_from_u64(42))
+    }
+}
+
+impl<T, R: Rng> KModesParams<T, R> {
+    /// Create a new K-Modes parameter builder with a custom RNG.
+    /// Defaults:
+    ///  * max_n_iterations: 100
+    ///  * n_runs: 10
+    ///  * init: KModesInit::Cao
+    ///  * verbose: false
+    pub fn new_with_rng(n_clusters: usize, rng: R) -> Self {
+        Self(KModesValidParams {
+            n_clusters,
+            max_n_iterations: 100,
+            n_runs: 10,
+            init: KModesInit::Cao,
+            verbose: false,
+            rng,
+        })
+    }
+
+    /// Set the maximum number of iterations for a single run.
+    pub fn max_n_iterations(mut self, max_n_iterations: u64) -> Self {
+        self.0.max_n_iterations = max_n_iterations;
+        self
+    }
+
+    /// Set the number of initialization runs (keeps the run with minimal cost/inertia).
+    pub fn n_runs(mut self, n_runs: usize) -> Self {
+        self.0.n_runs = n_runs;
+        self
+    }
+
+    /// Set the centroid initialization method (`Cao`, `Huang`, `Random`, or `Precomputed`).
+    pub fn init_method(mut self, init: KModesInit<T>) -> Self {
+        self.0.init = init;
+        self
+    }
+
+    /// Enable or disable verbose progress logging to stdout.
+    pub fn verbose(mut self, verbose: bool) -> Self {
+        self.0.verbose = verbose;
+        self
+    }
+
+    /// Set a custom random number generator.
+    pub fn with_rng<R2: Rng>(self, rng: R2) -> KModesParams<T, R2> {
+        KModesParams(KModesValidParams {
+            n_clusters: self.0.n_clusters,
+            max_n_iterations: self.0.max_n_iterations,
+            n_runs: self.0.n_runs,
+            init: self.0.init,
+            verbose: self.0.verbose,
+            rng,
+        })
+    }
+}
+
+impl<T, R: Rng> ParamGuard for KModesParams<T, R> {
+    type Checked = KModesValidParams<T, R>;
+    type Error = KModesParamsError;
+
+    fn check_ref(&self) -> Result<&Self::Checked, Self::Error> {
+        if self.0.n_clusters == 0 {
+            Err(KModesParamsError::NClusters)
+        } else if self.0.max_n_iterations == 0 {
+            Err(KModesParamsError::MaxIterations)
+        } else if self.0.n_runs == 0 {
+            Err(KModesParamsError::NRuns)
+        } else {
+            Ok(&self.0)
+        }
+    }
+
+    fn check(self) -> Result<Self::Checked, Self::Error> {
+        self.check_ref()?;
+        Ok(self.0)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use linfa::ParamGuard;
+    use ndarray_rand::rand::SeedableRng;
+    use rand_xoshiro::Xoshiro256Plus;
+
+    #[test]
+    fn autotraits() {
+        fn has_autotraits<T: Send + Sync + Sized + Unpin>() {}
+        has_autotraits::<KModesParamsError>();
+        has_autotraits::<KModesParams<usize, Xoshiro256Plus>>();
+        has_autotraits::<KModesValidParams<usize, Xoshiro256Plus>>();
+    }
+
+    #[test]
+    fn invalid_params() {
+        assert_eq!(
+            KModesParams::<usize, _>::new(0).check(),
+            Err(KModesParamsError::NClusters)
+        );
+        assert_eq!(
+            KModesParams::<usize, _>::new(2).max_n_iterations(0).check(),
+            Err(KModesParamsError::MaxIterations)
+        );
+        assert_eq!(
+            KModesParams::<usize, _>::new(2).n_runs(0).check(),
+            Err(KModesParamsError::NRuns)
+        );
+    }
+
+    #[test]
+    fn test_getters_and_builder_options() {
+        let params = KModesParams::<&str, _>::new(3)
+            .max_n_iterations(200)
+            .n_runs(15)
+            .init_method(KModesInit::Huang)
+            .verbose(true);
+
+        assert!(params.check_ref().is_ok());
+        let valid = params.check().unwrap();
+        assert_eq!(valid.n_clusters(), 3);
+        assert_eq!(valid.max_n_iterations(), 200);
+        assert_eq!(valid.n_runs(), 15);
+        assert_eq!(valid.init_method(), &KModesInit::Huang);
+        assert!(valid.verbose());
+        let _ = valid.rng();
+    }
+
+    #[test]
+    fn test_custom_rng() {
+        let rng1 = Xoshiro256Plus::seed_from_u64(12345);
+        let params1 = KModesParams::<usize, _>::new_with_rng(4, rng1);
+        let valid1 = params1.check().unwrap();
+        assert_eq!(valid1.n_clusters(), 4);
+
+        let rng2 = Xoshiro256Plus::seed_from_u64(67890);
+        let params2 = KModesParams::<usize, _>::new(2).with_rng(rng2);
+        let valid2 = params2.check().unwrap();
+        assert_eq!(valid2.n_clusters(), 2);
+    }
+
+    #[cfg(feature = "serde")]
+    #[cfg_attr(
+        feature = "serde",
+        derive(Serialize, Deserialize),
+        serde(crate = "serde_crate")
+    )]
+    #[derive(Clone, Debug, PartialEq)]
+    struct DummyRng;
+
+    #[cfg(feature = "serde")]
+    impl ndarray_rand::rand::RngCore for DummyRng {
+        fn next_u32(&mut self) -> u32 {
+            0
+        }
+        fn next_u64(&mut self) -> u64 {
+            0
+        }
+        fn fill_bytes(&mut self, dest: &mut [u8]) {
+            dest.fill(0);
+        }
+        fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), ndarray_rand::rand::Error> {
+            dest.fill(0);
+            Ok(())
+        }
+    }
+
+    #[cfg(feature = "serde")]
+    #[test]
+    fn test_serde_hyperparams() {
+        use ndarray_rand::rand::RngCore;
+        let params = KModesParams::<String, _>::new_with_rng(3, DummyRng)
+            .max_n_iterations(50)
+            .n_runs(5)
+            .init_method(KModesInit::Cao)
+            .verbose(false);
+        let valid = params.check().unwrap();
+
+        let serialized = serde_json::to_string(&valid).unwrap();
+        let deserialized: KModesValidParams<String, DummyRng> =
+            serde_json::from_str(&serialized).unwrap();
+
+        let mut dummy = DummyRng;
+        assert_eq!(dummy.next_u32(), 0);
+        assert_eq!(dummy.next_u64(), 0);
+        let mut buf = [0u8; 4];
+        dummy.fill_bytes(&mut buf);
+        assert!(dummy.try_fill_bytes(&mut buf).is_ok());
+
+        assert_eq!(valid.n_clusters(), deserialized.n_clusters());
+        assert_eq!(valid.max_n_iterations(), deserialized.max_n_iterations());
+        assert_eq!(valid.n_runs(), deserialized.n_runs());
+        assert_eq!(valid.init_method(), deserialized.init_method());
+        assert_eq!(valid.verbose(), deserialized.verbose());
+    }
+}
diff --git a/algorithms/linfa-clustering/src/k_modes/init.rs b/algorithms/linfa-clustering/src/k_modes/init.rs
new file mode 100644
index 000000000..a6226d9ef
--- /dev/null
+++ b/algorithms/linfa-clustering/src/k_modes/init.rs
@@ -0,0 +1,384 @@
+use ndarray::{Array2, ArrayView2};
+use ndarray_rand::rand::Rng;
+#[cfg(feature = "serde")]
+use serde_crate::{Deserialize, Serialize};
+use std::collections::{HashMap, HashSet};
+
+#[cfg_attr(
+    feature = "serde",
+    derive(Serialize, Deserialize),
+    serde(crate = "serde_crate")
+)]
+#[derive(Clone, Debug, PartialEq, Eq, Default)]
+#[non_exhaustive]
+/// Specifies the centroid initialization strategy for K-Modes.
+pub enum KModesInit<T> {
+    /// Density and dissimilarity-based initialization (Cao et al. \[2009\]). Default.
+    #[default]
+    Cao,
+    /// Density-based initialization (Huang \[1997, 1998\]).
+    Huang,
+    /// Randomly selects `n_clusters` unique observations from the dataset.
+    Random,
+    /// Precomputed initial centroids with shape `(n_clusters, n_features)`.
+    Precomputed(Array2<T>),
+}
+
+impl<T: crate::k_modes::EquivalenceTarget> KModesInit<T> {
+    /// Executes the configured initialization routine and returns initial centroids.
+    pub(crate) fn run<R: Rng>(
+        &self,
+        n_clusters: usize,
+        observations: ArrayView2<T>,
+        rng: &mut R,
+    ) -> Array2<T> {
+        match self {
+            Self::Cao => init_cao(observations, n_clusters),
+            Self::Huang => init_huang(observations, n_clusters, rng),
+            Self::Random => random_init(observations, n_clusters, rng),
+            Self::Precomputed(centroids) => {
+                assert_eq!(centroids.nrows(), n_clusters);
+                assert_eq!(centroids.ncols(), observations.ncols());
+                centroids.clone()
+            }
+        }
+    }
+}
+
+/// Cao initialization (Cao et al. \[2009\]): Density and dissimilarity-based selection.
+///
+/// Steps:
+/// 1. Compute marginal attribute density for every point $X_i$:
+///    $$\text{Dens}(X_i) = \sum_{j=1}^{M} \frac{\text{freq}_j(x_{ij})}{N \cdot M}$$
+///    where $\text{freq}_j(v)$ is the frequency of category $v$ across feature $j$,
+///    $N$ is the sample count, and $M$ is the feature count.
+/// 2. Pick the point with maximum density as the first centroid:
+///    $$C_0 = \text{argmax}_{i} \text{Dens}(X_i)$$
+/// 3. For subsequent centroids $k \in \{1, \dots, K-1\}$, choose the unselected point
+///    that maximizes the minimum density-weighted distance to existing centroids:
+///    $$C_k = \text{argmax}_{i \notin C} \left( \min_{c \in C} d(X_i, c) \cdot \text{Dens}(X_i) \right)$$
+///    This balances selecting high-density centers while maintaining spatial dispersion.
+pub(crate) fn init_cao<T: crate::k_modes::EquivalenceTarget>(
+    x: ArrayView2<T>,
+    n_clusters: usize,
+) -> Array2<T> {
+    let (n_points, n_attrs) = x.dim();
+    assert!(n_clusters <= n_points);
+
+    // Compute marginal attribute density per point
+    let mut dens = vec![0.0f64; n_points];
+    for iattr in 0..n_attrs {
+        let mut freq: HashMap<&T, usize> = HashMap::new();
+        for ipoint in 0..n_points {
+            *freq.entry(&x[[ipoint, iattr]]).or_insert(0) += 1;
+        }
+        let denom = (n_points * n_attrs) as f64;
+        for ipoint in 0..n_points {
+            let count = freq.get(&x[[ipoint, iattr]]).copied().unwrap_or(0);
+            dens[ipoint] += (count as f64) / denom;
+        }
+    }
+
+    let mut selected_indices = Vec::with_capacity(n_clusters);
+
+    // First centroid is the point with the highest density
+    let first_idx = dens
+        .iter()
+        .enumerate()
+        .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
+        .map(|(i, _)| i)
+        .unwrap_or(0);
+    selected_indices.push(first_idx);
+
+    // Subsequent centroids maximize the minimum density-weighted distance to chosen centroids
+    for _ in 1..n_clusters {
+        let mut best_idx = 0;
+        let mut best_score = -1.0f64;
+
+        for (ipoint, &point_dens) in dens.iter().enumerate() {
+            if selected_indices.contains(&ipoint) {
+                continue;
+            }
+
+            let row = x.row(ipoint);
+            let mut min_d = f64::INFINITY;
+
+            for &c_idx in &selected_indices {
+                let c_row = x.row(c_idx);
+                let dist = row.iter().zip(c_row.iter()).filter(|(a, b)| a != b).count() as f64;
+                let score = dist * point_dens;
+                if score < min_d {
+                    min_d = score;
+                }
+            }
+
+            if min_d > best_score {
+                best_score = min_d;
+                best_idx = ipoint;
+            }
+        }
+        selected_indices.push(best_idx);
+    }
+
+    Array2::from_shape_fn((n_clusters, n_attrs), |(i, j)| {
+        x[[selected_indices[i], j]].clone()
+    })
+}
+
+/// Huang initialization (Huang \[1997, 1998\]): Attribute frequency distribution sampling.
+///
+/// Steps:
+/// 1. Construct $K$ tentative synthetic centroids where each attribute $j$ is independently
+///    sampled proportional to its marginal category frequencies in the dataset.
+/// 2. Because synthetic combinations may not represent real observations or could cause empty clusters,
+///    each tentative centroid is mapped to its nearest distinct observation in $X$ via a linear scan.
+pub(crate) fn init_huang<T: crate::k_modes::EquivalenceTarget, R: Rng>(
+    x: ArrayView2<T>,
+    n_clusters: usize,
+    rng: &mut R,
+) -> Array2<T> {
+    let (n_points, n_attrs) = x.dim();
+    assert!(n_clusters <= n_points);
+
+    // Sample tentative centroids using column marginal distributions
+    let mut tentative = Array2::from_shape_fn((n_clusters, n_attrs), |(_, j)| {
+        let rand_idx = rng.gen_range(0..n_points);
+        x[[rand_idx, j]].clone()
+    });
+
+    // Replace each tentative centroid with the closest unselected observation in X
+    let mut selected_indices = HashSet::new();
+    for ik in 0..n_clusters {
+        let tent_row = tentative.row(ik);
+        let mut best_idx = 0;
+        let mut min_dist = usize::MAX;
+
+        for ipoint in 0..n_points {
+            if selected_indices.contains(&ipoint) {
+                continue;
+            }
+
+            let row = x.row(ipoint);
+            let dist = row
+                .iter()
+                .zip(tent_row.iter())
+                .filter(|(a, b)| a != b)
+                .count();
+
+            if dist < min_dist {
+                min_dist = dist;
+                best_idx = ipoint;
+            }
+        }
+
+        selected_indices.insert(best_idx);
+        for j in 0..n_attrs {
+            tentative[[ik, j]] = x[[best_idx, j]].clone();
+        }
+    }
+
+    tentative
+}
+
+/// Random initialization: selects `n_clusters` distinct data points without replacement.
+///
+/// Uses an ordered vector with a hash set to guarantee seed-deterministic output ordering.
+pub(crate) fn random_init<T: Clone, R: Rng>(
+    x: ArrayView2<T>,
+    n_clusters: usize,
+    rng: &mut R,
+) -> Array2<T> {
+    let (n_points, n_attrs) = x.dim();
+    assert!(n_clusters <= n_points);
+
+    let mut selected_indices = Vec::with_capacity(n_clusters);
+    let mut seen = HashSet::with_capacity(n_clusters);
+
+    while selected_indices.len() < n_clusters {
+        let idx = rng.gen_range(0..n_points);
+        if seen.insert(idx) {
+            selected_indices.push(idx);
+        }
+    }
+
+    Array2::from_shape_fn((n_clusters, n_attrs), |(i, j)| {
+        x[[selected_indices[i], j]].clone()
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ndarray::array;
+    use ndarray_rand::rand::SeedableRng;
+    use rand_xoshiro::Xoshiro256Plus;
+
+    #[test]
+    fn autotraits() {
+        fn has_autotraits<T: Send + Sync + Sized + Unpin>() {}
+        has_autotraits::<KModesInit<usize>>();
+    }
+
+    #[test]
+    fn test_init_cao_deterministic() {
+        let data = array![
+            ["A", "X", "1"],
+            ["A", "X", "1"],
+            ["A", "X", "2"],
+            ["B", "Y", "2"],
+            ["B", "Y", "2"],
+            ["C", "Z", "3"]
+        ];
+
+        let c1 = KModesInit::Cao.run(2, data.view(), &mut Xoshiro256Plus::seed_from_u64(42));
+        let c2 = KModesInit::Cao.run(2, data.view(), &mut Xoshiro256Plus::seed_from_u64(999));
+
+        assert_eq!(c1.dim(), (2, 3));
+        assert_eq!(c1, c2, "Cao initialization must be fully deterministic");
+
+        // The first centroid in Cao is the highest density point ("A", "X", "2")
+        let first_centroid = c1.row(0);
+        assert_eq!(first_centroid[0], "A");
+        assert_eq!(first_centroid[1], "X");
+        assert_eq!(first_centroid[2], "2");
+    }
+
+    #[test]
+    fn test_init_huang() {
+        let data = array![
+            ["A", "X"],
+            ["A", "X"],
+            ["B", "Y"],
+            ["B", "Y"],
+            ["C", "Z"],
+            ["C", "Z"]
+        ];
+
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        let centroids = KModesInit::Huang.run(2, data.view(), &mut rng);
+
+        assert_eq!(centroids.dim(), (2, 2));
+        // Verify centroids are drawn from valid data points
+        for row in centroids.rows() {
+            let matches = data.rows().into_iter().any(|d_row| d_row == row);
+            assert!(
+                matches,
+                "Huang centroid must match an existing data observation"
+            );
+        }
+    }
+
+    #[test]
+    fn test_init_random_unique_indices() {
+        let data = array![[1, 10], [2, 20], [3, 30], [4, 40], [5, 50]];
+
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        let centroids = KModesInit::Random.run(3, data.view(), &mut rng);
+
+        assert_eq!(centroids.dim(), (3, 2));
+        // Verify that 3 distinct rows from data were chosen
+        let r0 = centroids.row(0);
+        let r1 = centroids.row(1);
+        let r2 = centroids.row(2);
+        assert_ne!(r0, r1);
+        assert_ne!(r0, r2);
+        assert_ne!(r1, r2);
+    }
+
+    #[test]
+    fn test_init_precomputed() {
+        let data = array![["A", "X"], ["B", "Y"], ["C", "Z"]];
+        let precomputed = array![["A", "X"], ["C", "Z"]];
+
+        let init = KModesInit::Precomputed(precomputed.clone());
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        let centroids = init.run(2, data.view(), &mut rng);
+
+        assert_eq!(centroids, precomputed);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_init_precomputed_wrong_n_clusters() {
+        let data = array![["A", "X"], ["B", "Y"], ["C", "Z"]];
+        let precomputed = array![["A", "X"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        KModesInit::Precomputed(precomputed).run(2, data.view(), &mut rng);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_init_precomputed_wrong_features() {
+        let data = array![["A", "X"], ["B", "Y"], ["C", "Z"]];
+        let precomputed = array![["A", "X", "extra"], ["B", "Y", "extra"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        KModesInit::Precomputed(precomputed).run(2, data.view(), &mut rng);
+    }
+
+    #[test]
+    fn test_init_single_cluster() {
+        let data = array![["A", "X"], ["B", "Y"], ["A", "X"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+
+        let c_cao = KModesInit::Cao.run(1, data.view(), &mut rng);
+        assert_eq!(c_cao.dim(), (1, 2));
+
+        let c_huang = KModesInit::Huang.run(1, data.view(), &mut rng);
+        assert_eq!(c_huang.dim(), (1, 2));
+
+        let c_random = KModesInit::Random.run(1, data.view(), &mut rng);
+        assert_eq!(c_random.dim(), (1, 2));
+    }
+
+    #[test]
+    fn test_init_n_clusters_eq_n_points() {
+        let data = array![["A", "X"], ["B", "Y"], ["C", "Z"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+
+        let c_cao = KModesInit::Cao.run(3, data.view(), &mut rng);
+        assert_eq!(c_cao.dim(), (3, 2));
+
+        let c_huang = KModesInit::Huang.run(3, data.view(), &mut rng);
+        assert_eq!(c_huang.dim(), (3, 2));
+
+        let c_random = KModesInit::Random.run(3, data.view(), &mut rng);
+        assert_eq!(c_random.dim(), (3, 2));
+    }
+
+    #[cfg(feature = "serde")]
+    #[test]
+    fn test_serde_init() {
+        let init: KModesInit<String> = KModesInit::Cao;
+        let serialized = serde_json::to_string(&init).unwrap();
+        let deserialized: KModesInit<String> = serde_json::from_str(&serialized).unwrap();
+        assert_eq!(init, deserialized);
+
+        let precomputed = KModesInit::Precomputed(array![["A".to_string(), "B".to_string()]]);
+        let serialized_pre = serde_json::to_string(&precomputed).unwrap();
+        let deserialized_pre: KModesInit<String> = serde_json::from_str(&serialized_pre).unwrap();
+        assert_eq!(precomputed, deserialized_pre);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_init_cao_n_clusters_exceeds_points() {
+        let data = array![["A", "1"]];
+        init_cao(data.view(), 5);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_init_huang_n_clusters_exceeds_points() {
+        let data = array![["A", "1"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        init_huang(data.view(), 5, &mut rng);
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_init_random_n_clusters_exceeds_points() {
+        let data = array![["A", "1"]];
+        let mut rng = Xoshiro256Plus::seed_from_u64(42);
+        random_init(data.view(), 5, &mut rng);
+    }
+}
diff --git a/algorithms/linfa-clustering/src/k_modes/mod.rs b/algorithms/linfa-clustering/src/k_modes/mod.rs
new file mode 100644
index 000000000..f9616d989
--- /dev/null
+++ b/algorithms/linfa-clustering/src/k_modes/mod.rs
@@ -0,0 +1,9 @@
+mod algorithm;
+mod errors;
+mod hyperparams;
+mod init;
+
+pub use algorithm::*;
+pub use errors::*;
+pub use hyperparams::*;
+pub use init::*;
diff --git a/algorithms/linfa-clustering/src/lib.rs b/algorithms/linfa-clustering/src/lib.rs
index 7418bc64a..cab61b534 100644
--- a/algorithms/linfa-clustering/src/lib.rs
+++ b/algorithms/linfa-clustering/src/lib.rs
@@ -14,6 +14,7 @@
 //!
 //! Right now `linfa-clustering` provides the following clustering algorithms:
 //! * [K-Means](KMeans)
+//! * [K-Modes](KModes)
 //! * [DBSCAN](Dbscan)
 //! * [Approximated DBSCAN](AppxDbscan) (Currently an alias for DBSCAN, due to its superior
 //!   performance)
@@ -25,11 +26,13 @@ mod dbscan;
 mod gaussian_mixture;
 #[allow(clippy::new_ret_no_self)]
 mod k_means;
+mod k_modes;
 mod optics;
 
 pub use dbscan::*;
 pub use gaussian_mixture::*;
 pub use k_means::*;
+pub use k_modes::*;
 pub use optics::*;
 
 // Approx DBSCAN is currently an alias for DBSCAN, due to the old Approx DBSCAN implementation's