From 1457dca42eb881562e467e0af6fe9a796b1ebd82 Mon Sep 17 00:00:00 2001 From: Ricardo Carvalho Santos Date: Sat, 11 Jul 2026 09:07:08 -0300 Subject: [PATCH 1/5] Add Kernel PCA decomposition Implements Scholar.Decomposition.KernelPCA, the kernel extension of PCA that runs the decomposition in a reproducing kernel Hilbert space, capturing non-linear structure that ordinary PCA cannot. Supports the linear, poly, rbf, sigmoid and cosine kernels, with fit/2, transform/2 and fit_transform/2 mirroring the existing PCA API. The kernel matrix is double-centered, decomposed with Nx.LinAlg.eigh, and the eigenvectors are sign-flipped for deterministic output. Results match scikit-learn's KernelPCA on all five kernels (covered by the tests). Addresses the Kernel PCA item of #246. --- lib/scholar/decomposition/kernel_pca.ex | 268 ++++++++++++++++++ mix.exs | 1 + .../scholar/decomposition/kernel_pca_test.exs | 194 +++++++++++++ 3 files changed, 463 insertions(+) create mode 100644 lib/scholar/decomposition/kernel_pca.ex create mode 100644 test/scholar/decomposition/kernel_pca_test.exs diff --git a/lib/scholar/decomposition/kernel_pca.ex b/lib/scholar/decomposition/kernel_pca.ex new file mode 100644 index 00000000..679f9dca --- /dev/null +++ b/lib/scholar/decomposition/kernel_pca.ex @@ -0,0 +1,268 @@ +defmodule Scholar.Decomposition.KernelPCA do + @moduledoc """ + Kernel Principal Component Analysis (Kernel PCA). + + Kernel PCA is an extension of `Scholar.Decomposition.PCA` that performs the + decomposition in a reproducing kernel Hilbert space instead of the original + feature space. Replacing the inner product with a kernel function lets the + method capture non-linear structure that ordinary PCA cannot. + + The time complexity is $O(N^3)$ where $N$ is the number of samples, since it + relies on the eigendecomposition of the $N \\times N$ kernel matrix. + + References: + + * [1] Schölkopf, B., Smola, A., & Müller, K. R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5), 1299-1319. + """ + import Nx.Defn + + @derive {Nx.Container, + keep: [:kernel, :gamma, :degree, :coef0], + containers: [:eigenvalues, :eigenvectors, :x_fit, :kernel_fit_rows, :kernel_fit_all]} + defstruct [ + :eigenvalues, + :eigenvectors, + :x_fit, + :kernel_fit_rows, + :kernel_fit_all, + :kernel, + :gamma, + :degree, + :coef0 + ] + + opts = [ + num_components: [ + required: true, + type: :pos_integer, + doc: "The number of principal components to keep." + ], + kernel: [ + type: {:in, [:linear, :poly, :rbf, :sigmoid, :cosine]}, + default: :linear, + doc: "The kernel used to compute the pairwise similarities between samples." + ], + gamma: [ + type: {:or, [:float, nil]}, + default: nil, + doc: """ + Kernel coefficient for the `:rbf`, `:poly` and `:sigmoid` kernels. + When `nil` it defaults to `1 / num_features`. + """ + ], + degree: [ + type: :pos_integer, + default: 3, + doc: "Degree of the `:poly` kernel." + ], + coef0: [ + type: :float, + default: 1.0, + doc: "Independent term of the `:poly` and `:sigmoid` kernels." + ] + ] + + @opts_schema NimbleOptions.new!(opts) + + @doc """ + Fits a Kernel PCA for sample inputs `x`. + + ## Options + + #{NimbleOptions.docs(@opts_schema)} + + ## Return Values + + The function returns a struct with the following parameters: + + * `:eigenvalues` - The eigenvalues of the centered kernel matrix, sorted in + decreasing order. + + * `:eigenvectors` - The eigenvectors of the centered kernel matrix. + + * `:x_fit` - The training data, kept to compute the kernel against new samples. + + * `:kernel_fit_rows` - Column means of the training kernel matrix, used to + center the kernel of new samples. + + * `:kernel_fit_all` - Mean of the whole training kernel matrix. + + ## Examples + + iex> x = Nx.tensor([[0.5, 0.2, 0.8], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) + iex> kpca = Scholar.Decomposition.KernelPCA.fit(x, num_components: 2, kernel: :rbf) + iex> kpca.eigenvalues + Nx.tensor([0.3644832372665405, 0.26074573397636414]) + """ + deftransform fit(x, opts \\ []) do + opts = NimbleOptions.validate!(opts, @opts_schema) + + if Nx.rank(x) != 2 do + raise ArgumentError, + """ + expected input tensor to have shape {num_samples, num_features}, \ + got tensor with shape: #{inspect(Nx.shape(x))}\ + """ + end + + {num_samples, num_features} = Nx.shape(x) + num_components = opts[:num_components] + + if num_components > num_samples do + raise ArgumentError, + """ + num_components must be less than or equal to \ + num_samples = #{num_samples}, got #{num_components}\ + """ + end + + opts = Keyword.put(opts, :gamma, opts[:gamma] || 1.0 / num_features) + fit_n(x, opts) + end + + defnp fit_n(x, opts) do + kernel = kernel_matrix(x, x, opts) + {kernel_centered, kernel_fit_rows, kernel_fit_all} = center_fit(kernel) + + {eigenvalues, eigenvectors} = top_eigen(kernel_centered, opts[:num_components]) + + %__MODULE__{ + eigenvalues: eigenvalues, + eigenvectors: eigenvectors, + x_fit: x, + kernel_fit_rows: kernel_fit_rows, + kernel_fit_all: kernel_fit_all, + kernel: opts[:kernel], + gamma: opts[:gamma], + degree: opts[:degree], + coef0: opts[:coef0] + } + end + + @doc """ + For a fitted `model` projects samples `x` onto the principal components. + + ## Return Values + + The function returns a tensor with the decomposed data. + + ## Examples + + iex> x = Nx.tensor([[0.5, 0.2, 0.8], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) + iex> kpca = Scholar.Decomposition.KernelPCA.fit(x, num_components: 2, kernel: :rbf) + iex> Scholar.Decomposition.KernelPCA.transform(kpca, Nx.tensor([[0.5, 0.5, 0.5]])) + Nx.tensor( + [ + [0.12500189244747162, 0.029097510501742363] + ] + ) + """ + deftransform transform(%__MODULE__{} = model, x) do + if Nx.rank(x) != 2 do + raise ArgumentError, + """ + expected input tensor to have shape {num_samples, num_features}, \ + got tensor with shape: #{inspect(Nx.shape(x))}\ + """ + end + + transform_n(model, x) + end + + defnp transform_n(%__MODULE__{x_fit: x_fit} = model, x) do + kernel = kernel_matrix(x, x_fit, model_opts(model)) + kernel_centered = center_transform(kernel, model.kernel_fit_rows, model.kernel_fit_all) + project(kernel_centered, model.eigenvalues, model.eigenvectors) + end + + @doc """ + Fits a Kernel PCA on `x` and projects `x` onto the principal components. + + ## Options + + #{NimbleOptions.docs(@opts_schema)} + + ## Return Values + + The function returns a tensor with the decomposed data. + + ## Examples + + iex> x = Nx.tensor([[0.5, 0.2, 0.8], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) + iex> Scholar.Decomposition.KernelPCA.fit_transform(x, num_components: 2, kernel: :rbf) + Nx.tensor( + [ + [-0.13561560213565826, -0.16519638895988464], + [0.021600963547825813, 0.44114428758621216], + [0.4687873125076294, -0.15764620900154114], + [-0.35477250814437866, -0.11830171197652817] + ] + ) + """ + deftransform fit_transform(x, opts \\ []) do + fit(x, opts) |> project_fit() + end + + defnp project_fit(model) do + model.eigenvectors * Nx.sqrt(model.eigenvalues) + end + + # Kernel matrix between `x` (m samples) and `y` (n samples), shape {m, n}. + defnp kernel_matrix(x, y, opts) do + case opts[:kernel] do + :linear -> + Nx.dot(x, [1], y, [1]) + + :poly -> + (opts[:gamma] * Nx.dot(x, [1], y, [1]) + opts[:coef0]) ** opts[:degree] + + :rbf -> + Nx.exp(-opts[:gamma] * Scholar.Metrics.Distance.pairwise_squared_euclidean(x, y)) + + :sigmoid -> + Nx.tanh(opts[:gamma] * Nx.dot(x, [1], y, [1]) + opts[:coef0]) + + :cosine -> + x_normalized = x / Nx.sqrt(Nx.sum(x * x, axes: [1], keep_axes: true)) + y_normalized = y / Nx.sqrt(Nx.sum(y * y, axes: [1], keep_axes: true)) + Nx.dot(x_normalized, [1], y_normalized, [1]) + end + end + + deftransformp model_opts(model) do + [kernel: model.kernel, gamma: model.gamma, degree: model.degree, coef0: model.coef0] + end + + # Double centering of the training kernel, plus the statistics reused at transform time. + defnp center_fit(kernel) do + column_means = Nx.mean(kernel, axes: [0], keep_axes: true) + row_means = Nx.mean(kernel, axes: [1], keep_axes: true) + all_mean = Nx.mean(kernel) + centered = kernel - column_means - row_means + all_mean + {centered, Nx.squeeze(column_means, axes: [0]), all_mean} + end + + defnp center_transform(kernel, kernel_fit_rows, kernel_fit_all) do + predict_row_means = Nx.mean(kernel, axes: [1], keep_axes: true) + kernel - kernel_fit_rows - predict_row_means + kernel_fit_all + end + + # Eigenvectors of the centered kernel, sorted by decreasing eigenvalue and + # sign-flipped so the largest absolute entry of each vector is positive. + defnp top_eigen(kernel_centered, num_components) do + {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(kernel_centered, eps: 1.0e-8) + + order = Nx.argsort(eigenvalues, direction: :desc) + eigenvalues = Nx.take(eigenvalues, order)[0..(num_components - 1)] + eigenvectors = Nx.take(eigenvectors, order, axis: 1)[[.., 0..(num_components - 1)]] + + max_abs = eigenvectors |> Nx.abs() |> Nx.argmax(axis: 0, keep_axis: true) + signs = eigenvectors |> Nx.take_along_axis(max_abs, axis: 0) |> Nx.sign() + {eigenvalues, eigenvectors * signs} + end + + defnp project(kernel_centered, eigenvalues, eigenvectors) do + scaled_eigenvectors = eigenvectors / Nx.sqrt(eigenvalues) + Nx.dot(kernel_centered, [1], scaled_eigenvectors, [0]) + end +end diff --git a/mix.exs b/mix.exs index c57c855c..9428cb46 100644 --- a/mix.exs +++ b/mix.exs @@ -72,6 +72,7 @@ defmodule Scholar.MixProject do Scholar.Cluster.GaussianMixture, Scholar.Cluster.Hierarchical, Scholar.Cluster.KMeans, + Scholar.Decomposition.KernelPCA, Scholar.Decomposition.PCA, Scholar.Integrate, Scholar.Interpolation.BezierSpline, diff --git a/test/scholar/decomposition/kernel_pca_test.exs b/test/scholar/decomposition/kernel_pca_test.exs new file mode 100644 index 00000000..3000cd8c --- /dev/null +++ b/test/scholar/decomposition/kernel_pca_test.exs @@ -0,0 +1,194 @@ +defmodule Scholar.Decomposition.KernelPCATest do + use Scholar.Case, async: true + alias Scholar.Decomposition.KernelPCA + doctest KernelPCA + + defp x do + Nx.tensor([ + [0.5, 0.2, 0.8], + [1.0, 0.5, 0.2], + [0.3, 1.0, 0.7], + [0.9, 0.1, 1.0], + [0.4, 0.8, 0.3], + [0.6, 0.2, 0.9] + ]) + end + + defp x_test do + Nx.tensor([[0.5, 0.5, 0.5], [0.1, 0.9, 0.2]]) + end + + # Reference values taken from scikit-learn (sklearn.decomposition.KernelPCA). + describe "linear kernel" do + test "fit/2 eigenvalues" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :linear) + assert_all_close(model.eigenvalues, Nx.tensor([1.02576507, 0.49150841])) + end + + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :linear) + + assert_all_close( + z, + Nx.tensor([ + [0.2480282, -0.14309686], + [-0.13098757, 0.57609125], + [-0.49936002, -0.3421911], + [0.56523378, 0.01348409], + [-0.5166106, 0.03850971], + [0.33369622, -0.14279708] + ]), + atol: 1.0e-3 + ) + end + + test "transform/2" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :linear) + + assert_all_close( + KernelPCA.transform(model, x_test()), + Nx.tensor([[-0.1434428, 0.01749365], [-0.74814817, -0.11777146]]), + atol: 1.0e-3 + ) + end + end + + describe "rbf kernel" do + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :rbf) + + assert_all_close(model_eigenvalues(:rbf), Nx.tensor([0.57052273, 0.28089124])) + + assert_all_close( + z, + Nx.tensor([ + [0.19920416, -0.10217718], + [-0.11443324, 0.43767286], + [-0.3679805, -0.26212427], + [0.40563184, 0.01272026], + [-0.38548342, 0.01315336], + [0.26306115, -0.09924505] + ]), + atol: 1.0e-3 + ) + end + + test "transform/2" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :rbf) + + assert_all_close( + KernelPCA.transform(model, x_test()), + Nx.tensor([[-0.11224132, 0.01060868], [-0.49994934, -0.10212445]]), + atol: 1.0e-3 + ) + end + end + + describe "poly kernel" do + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :poly, degree: 3, coef0: 1.0) + + assert_all_close( + z, + Nx.tensor([ + [0.27508957, -0.15923293], + [-0.16730544, 0.7908652], + [-0.73099016, -0.47917877], + [0.85131884, -0.04410387], + [-0.65531824, 0.08343613], + [0.42720543, -0.19178576] + ]), + atol: 1.0e-3 + ) + end + end + + describe "sigmoid kernel" do + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :sigmoid, coef0: 1.0) + + assert_all_close( + z, + Nx.tensor([ + [-0.08427136, -0.04391895], + [0.04163971, 0.15793664], + [0.12574582, -0.09492698], + [-0.13989375, 0.01645784], + [0.15586995, 0.00172166], + [-0.09909037, -0.03727021] + ]), + atol: 1.0e-3 + ) + end + end + + describe "cosine kernel" do + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :cosine) + + assert_all_close( + z, + Nx.tensor([ + [-0.30958794, -0.13265305], + [0.1569742, 0.49544337], + [0.37577318, -0.30154498], + [-0.40180281, 0.06364262], + [0.50860284, -0.02523757], + [-0.32995947, -0.09965038] + ]), + atol: 1.0e-3 + ) + end + end + + describe "general behaviour" do + test "fit_transform/2 equals fit/2 followed by transform/2 on training data" do + for kernel <- [:linear, :poly, :rbf, :sigmoid, :cosine] do + opts = [num_components: 2, kernel: kernel] + model = KernelPCA.fit(x(), opts) + + assert_all_close( + KernelPCA.fit_transform(x(), opts), + KernelPCA.transform(model, x()), + atol: 1.0e-3 + ) + end + end + + test "the number of components sets the output dimension" do + z = KernelPCA.fit_transform(x(), num_components: 3, kernel: :rbf) + assert Nx.shape(z) == {6, 3} + end + + test "fit/2 and transform/2 work with jit_apply" do + model = Nx.Defn.jit_apply(&KernelPCA.fit(&1, num_components: 2, kernel: :rbf), [x()]) + z = Nx.Defn.jit_apply(&KernelPCA.transform/2, [model, x_test()]) + assert Nx.shape(z) == {2, 2} + end + + test "fit/2 propagates input precision (f64)" do + x = Nx.as_type(x(), :f64) + model = KernelPCA.fit(x, num_components: 2, kernel: :rbf) + assert Nx.type(model.eigenvalues) == {:f, 64} + assert Nx.type(KernelPCA.transform(model, Nx.as_type(x_test(), :f64))) == {:f, 64} + end + end + + describe "input validation" do + test "fit/2 requires a rank-2 tensor" do + assert_raise ArgumentError, + ~r/expected input tensor to have shape \{num_samples, num_features\}/, + fn -> KernelPCA.fit(Nx.iota({5}), num_components: 2) end + end + + test "num_components may not exceed num_samples" do + assert_raise ArgumentError, + "num_components must be less than or equal to num_samples = 6, got 7", + fn -> KernelPCA.fit(x(), num_components: 7) end + end + end + + defp model_eigenvalues(kernel) do + KernelPCA.fit(x(), num_components: 2, kernel: kernel).eigenvalues + end +end From 0bc738b7ea7334ff645098cf4c49936009ec5dc7 Mon Sep 17 00:00:00 2001 From: Ricardo Carvalho Santos Date: Sat, 11 Jul 2026 13:38:42 -0300 Subject: [PATCH 2/5] Fix robustness issues found in Kernel PCA review Three issues found while reviewing the initial implementation: * transform/2 now validates the input has the same number of features used to fit the model, matching PCA's own validation, instead of raising an Nx-internal shape error. * The cosine kernel no longer divides by zero for all-zero rows; it treats them as zero similarity, matching scikit-learn's normalize(). * Centering the kernel matrix can leave it asymmetric by a few ULPs, which Nx.LinAlg.eigh rejects outright since it requires exact symmetry (this could crash fit/2 depending on kernel/gamma/dataset). The centered kernel is now explicitly symmetrized before the decomposition. Also expands test coverage: transform/2 and eigenvalues for poly/sigmoid/ cosine (previously only checked via fit_transform/2), custom gamma and degree, the num_components == num_samples boundary, default kernel/gamma, eigenvalue ordering, and regression tests for the three fixes above. --- lib/scholar/decomposition/kernel_pca.ex | 28 +++- .../scholar/decomposition/kernel_pca_test.exs | 145 ++++++++++++++++++ 2 files changed, 170 insertions(+), 3 deletions(-) diff --git a/lib/scholar/decomposition/kernel_pca.ex b/lib/scholar/decomposition/kernel_pca.ex index 679f9dca..f53ea53f 100644 --- a/lib/scholar/decomposition/kernel_pca.ex +++ b/lib/scholar/decomposition/kernel_pca.ex @@ -166,6 +166,19 @@ defmodule Scholar.Decomposition.KernelPCA do """ end + num_features_seen = Nx.axis_size(model.x_fit, 1) + num_features = Nx.axis_size(x, 1) + + if num_features_seen != num_features do + raise ArgumentError, + """ + expected input tensor to have the same number of features \ + as tensor used to fit the model, \ + got #{inspect(num_features)} \ + and #{inspect(num_features_seen)}\ + """ + end + transform_n(model, x) end @@ -223,12 +236,18 @@ defmodule Scholar.Decomposition.KernelPCA do Nx.tanh(opts[:gamma] * Nx.dot(x, [1], y, [1]) + opts[:coef0]) :cosine -> - x_normalized = x / Nx.sqrt(Nx.sum(x * x, axes: [1], keep_axes: true)) - y_normalized = y / Nx.sqrt(Nx.sum(y * y, axes: [1], keep_axes: true)) + # rows with zero norm are left as all-zeros instead of dividing by zero + x_normalized = x / safe_norm(x) + y_normalized = y / safe_norm(y) Nx.dot(x_normalized, [1], y_normalized, [1]) end end + defnp safe_norm(x) do + norm = Nx.sqrt(Nx.sum(x * x, axes: [1], keep_axes: true)) + Nx.select(norm == 0, 1.0, norm) + end + deftransformp model_opts(model) do [kernel: model.kernel, gamma: model.gamma, degree: model.degree, coef0: model.coef0] end @@ -250,7 +269,10 @@ defmodule Scholar.Decomposition.KernelPCA do # Eigenvectors of the centered kernel, sorted by decreasing eigenvalue and # sign-flipped so the largest absolute entry of each vector is positive. defnp top_eigen(kernel_centered, num_components) do - {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(kernel_centered, eps: 1.0e-8) + # Centering can leave tiny floating-point asymmetries, but eigh requires + # exact symmetry, so it is enforced explicitly before the decomposition. + symmetric_kernel = (kernel_centered + Nx.transpose(kernel_centered)) / 2 + {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(symmetric_kernel, eps: 1.0e-8) order = Nx.argsort(eigenvalues, direction: :desc) eigenvalues = Nx.take(eigenvalues, order)[0..(num_components - 1)] diff --git a/test/scholar/decomposition/kernel_pca_test.exs b/test/scholar/decomposition/kernel_pca_test.exs index 3000cd8c..ec1035e6 100644 --- a/test/scholar/decomposition/kernel_pca_test.exs +++ b/test/scholar/decomposition/kernel_pca_test.exs @@ -101,6 +101,38 @@ defmodule Scholar.Decomposition.KernelPCATest do atol: 1.0e-3 ) end + + test "fit/2 eigenvalues" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :poly, degree: 3, coef0: 1.0) + assert_all_close(model.eigenvalues, Nx.tensor([1.97470225, 0.9261237])) + end + + test "transform/2" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :poly, degree: 3, coef0: 1.0) + + assert_all_close( + KernelPCA.transform(model, x_test()), + Nx.tensor([[-0.21107767, 0.05450511], [-0.87419523, -0.07912082]]), + atol: 1.0e-3 + ) + end + + test "custom degree" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :poly, degree: 2, coef0: 1.0) + + assert_all_close( + z, + Nx.tensor([ + [0.21428251, -0.12551537], + [-0.11978673, 0.55138703], + [-0.49389717, -0.32901762], + [0.56577048, -0.00833521], + [-0.47484395, 0.04791742], + [0.30847487, -0.13643626] + ]), + atol: 1.0e-3 + ) + end end describe "sigmoid kernel" do @@ -120,6 +152,21 @@ defmodule Scholar.Decomposition.KernelPCATest do atol: 1.0e-3 ) end + + test "fit/2 eigenvalues" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :sigmoid, coef0: 1.0) + assert_all_close(model.eigenvalues, Nx.tensor([0.07833214, 0.03754688])) + end + + test "transform/2" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :sigmoid, coef0: 1.0) + + assert_all_close( + KernelPCA.transform(model, x_test()), + Nx.tensor([[0.03661575, -0.00068769], [0.24706743, -0.05864022]]), + atol: 1.0e-3 + ) + end end describe "cosine kernel" do @@ -139,6 +186,67 @@ defmodule Scholar.Decomposition.KernelPCATest do atol: 1.0e-3 ) end + + test "fit/2 eigenvalues" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :cosine) + assert_all_close(model.eigenvalues, Nx.tensor([0.79068667, 0.36860785])) + end + + test "transform/2" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :cosine) + + assert_all_close( + KernelPCA.transform(model, x_test()), + Nx.tensor([[0.11596968, -0.0139805], [0.72969734, -0.22251837]]), + atol: 1.0e-3 + ) + end + + test "rows with zero norm do not produce NaN" do + x_with_zero_row = + Nx.tensor([[0.0, 0.0, 0.0], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) + + model = KernelPCA.fit(x_with_zero_row, num_components: 2, kernel: :cosine) + + # Reference values taken from scikit-learn, which treats zero-norm rows + # as an all-zero vector instead of dividing by zero. + assert_all_close(model.eigenvalues, Nx.tensor([0.59310028, 0.38942085])) + + assert_all_close( + KernelPCA.fit_transform(x_with_zero_row, num_components: 2, kernel: :cosine), + Nx.tensor([ + [0.66198055, -0.06451232], + [-0.2645941, -0.1783384], + [-0.14431492, 0.52389659], + [-0.25307153, -0.28104588] + ]), + atol: 1.0e-3 + ) + end + end + + describe "rbf kernel with custom gamma" do + test "fit_transform/2" do + z = KernelPCA.fit_transform(x(), num_components: 2, kernel: :rbf, gamma: 2.0) + + assert_all_close( + z, + Nx.tensor([ + [-0.43668454, -0.12918044], + [0.29479921, 0.79022216], + [0.5971787, -0.46764162], + [-0.56175187, 0.01788282], + [0.63960914, -0.09610005], + [-0.53315064, -0.11518287] + ]), + atol: 1.0e-3 + ) + end + + test "fit/2 keeps the given gamma instead of the default" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :rbf, gamma: 2.0) + assert model.gamma == 2.0 + end end describe "general behaviour" do @@ -172,6 +280,26 @@ defmodule Scholar.Decomposition.KernelPCATest do assert Nx.type(model.eigenvalues) == {:f, 64} assert Nx.type(KernelPCA.transform(model, Nx.as_type(x_test(), :f64))) == {:f, 64} end + + test "kernel defaults to :linear" do + assert KernelPCA.fit(x(), num_components: 2).kernel == :linear + end + + test "gamma defaults to 1 / num_features" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :rbf) + assert_all_close(Nx.tensor(model.gamma), Nx.tensor(1.0 / 3)) + end + + test "num_components may equal num_samples" do + model = KernelPCA.fit(x(), num_components: 6, kernel: :linear) + assert Nx.shape(model.eigenvalues) == {6} + end + + test "eigenvalues are sorted in decreasing order" do + model = KernelPCA.fit(x(), num_components: 4, kernel: :rbf) + sorted = model.eigenvalues |> Nx.to_flat_list() |> Enum.sort(:desc) + assert Nx.to_flat_list(model.eigenvalues) == sorted + end end describe "input validation" do @@ -186,6 +314,23 @@ defmodule Scholar.Decomposition.KernelPCATest do "num_components must be less than or equal to num_samples = 6, got 7", fn -> KernelPCA.fit(x(), num_components: 7) end end + + test "transform/2 requires a rank-2 tensor" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :linear) + + assert_raise ArgumentError, + ~r/expected input tensor to have shape \{num_samples, num_features\}/, + fn -> KernelPCA.transform(model, Nx.iota({3})) end + end + + test "transform/2 requires the same number of features used to fit the model" do + model = KernelPCA.fit(x(), num_components: 2, kernel: :linear) + + assert_raise ArgumentError, + "expected input tensor to have the same number of features " <> + "as tensor used to fit the model, got 2 and 3", + fn -> KernelPCA.transform(model, Nx.tensor([[0.1, 0.2]])) end + end end defp model_eigenvalues(kernel) do From 95360cf82a072daad818c01c3a3a9ca66fe3a76e Mon Sep 17 00:00:00 2001 From: Ricardo Carvalho Santos Date: Sat, 11 Jul 2026 20:20:14 -0300 Subject: [PATCH 3/5] Handle zero eigenvalues and match eigh tolerance to input precision Verifying every output against scikit-learn surfaced two more issues: * With num_components close to num_samples, the smallest eigenvalue of the centered kernel is zero up to floating-point noise (the centering forces it), and could come out slightly negative, crashing Nx.sqrt in fit_transform/2. Indefinite kernels such as sigmoid can also produce genuinely negative eigenvalues. They are now clipped to zero and their components project to zero, as in scikit-learn. * The eps previously passed to Nx.LinAlg.eigh (1.0e-8) stopped the iteration too early for f64 inputs, leaving eigenpair residuals around 1.0e-4 where LAPACK reaches 1.0e-16, while tighter tolerances made f32 accumulate rounding noise past convergence. The tolerance now matches the input precision (1.0e-8 for f32, 1.0e-11 for f64), which brings the f64 projections within 1.0e-5 of the SciPy reference in the worst-conditioned case measured (500x better than before). Also validates gamma as a positive number and accepts integers for gamma and coef0, following the option style used elsewhere in Scholar. --- lib/scholar/decomposition/kernel_pca.ex | 27 ++++++++++++++----- .../scholar/decomposition/kernel_pca_test.exs | 27 +++++++++++++++++-- 2 files changed, 46 insertions(+), 8 deletions(-) diff --git a/lib/scholar/decomposition/kernel_pca.ex b/lib/scholar/decomposition/kernel_pca.ex index f53ea53f..ab104b41 100644 --- a/lib/scholar/decomposition/kernel_pca.ex +++ b/lib/scholar/decomposition/kernel_pca.ex @@ -43,7 +43,7 @@ defmodule Scholar.Decomposition.KernelPCA do doc: "The kernel used to compute the pairwise similarities between samples." ], gamma: [ - type: {:or, [:float, nil]}, + type: {:or, [{:custom, Scholar.Options, :positive_number, []}, nil]}, default: nil, doc: """ Kernel coefficient for the `:rbf`, `:poly` and `:sigmoid` kernels. @@ -56,7 +56,7 @@ defmodule Scholar.Decomposition.KernelPCA do doc: "Degree of the `:poly` kernel." ], coef0: [ - type: :float, + type: {:or, [:float, :integer]}, default: 1.0, doc: "Independent term of the `:poly` and `:sigmoid` kernels." ] @@ -117,6 +117,13 @@ defmodule Scholar.Decomposition.KernelPCA do end opts = Keyword.put(opts, :gamma, opts[:gamma] || 1.0 / num_features) + + # The stopping criterion of Nx.LinAlg.eigh must match the input precision: + # a tolerance below the floating-point resolution makes the iteration run + # long past convergence and accumulate rounding noise instead. + eigh_eps = if Nx.type(x) == {:f, 64}, do: 1.0e-11, else: 1.0e-8 + opts = Keyword.put(opts, :eigh_eps, eigh_eps) + fit_n(x, opts) end @@ -124,7 +131,8 @@ defmodule Scholar.Decomposition.KernelPCA do kernel = kernel_matrix(x, x, opts) {kernel_centered, kernel_fit_rows, kernel_fit_all} = center_fit(kernel) - {eigenvalues, eigenvectors} = top_eigen(kernel_centered, opts[:num_components]) + {eigenvalues, eigenvectors} = + top_eigen(kernel_centered, opts[:num_components], opts[:eigh_eps]) %__MODULE__{ eigenvalues: eigenvalues, @@ -268,23 +276,30 @@ defmodule Scholar.Decomposition.KernelPCA do # Eigenvectors of the centered kernel, sorted by decreasing eigenvalue and # sign-flipped so the largest absolute entry of each vector is positive. - defnp top_eigen(kernel_centered, num_components) do + defnp top_eigen(kernel_centered, num_components, eigh_eps) do # Centering can leave tiny floating-point asymmetries, but eigh requires # exact symmetry, so it is enforced explicitly before the decomposition. symmetric_kernel = (kernel_centered + Nx.transpose(kernel_centered)) / 2 - {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(symmetric_kernel, eps: 1.0e-8) + {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(symmetric_kernel, eps: eigh_eps) order = Nx.argsort(eigenvalues, direction: :desc) eigenvalues = Nx.take(eigenvalues, order)[0..(num_components - 1)] eigenvectors = Nx.take(eigenvectors, order, axis: 1)[[.., 0..(num_components - 1)]] + # negative eigenvalues carry no usable projection, either floating-point + # noise around zero (the centering itself forces one zero eigenvalue) or + # produced by an indefinite kernel such as sigmoid, so they are clipped + eigenvalues = Nx.max(eigenvalues, 0) + max_abs = eigenvectors |> Nx.abs() |> Nx.argmax(axis: 0, keep_axis: true) signs = eigenvectors |> Nx.take_along_axis(max_abs, axis: 0) |> Nx.sign() {eigenvalues, eigenvectors * signs} end + # components with a zero eigenvalue are projected to zero, as in scikit-learn defnp project(kernel_centered, eigenvalues, eigenvectors) do - scaled_eigenvectors = eigenvectors / Nx.sqrt(eigenvalues) + safe_eigenvalues = Nx.select(eigenvalues == 0, 1.0, eigenvalues) + scaled_eigenvectors = eigenvectors / Nx.sqrt(safe_eigenvalues) * (eigenvalues > 0) Nx.dot(kernel_centered, [1], scaled_eigenvectors, [0]) end end diff --git a/test/scholar/decomposition/kernel_pca_test.exs b/test/scholar/decomposition/kernel_pca_test.exs index ec1035e6..15f1e9e7 100644 --- a/test/scholar/decomposition/kernel_pca_test.exs +++ b/test/scholar/decomposition/kernel_pca_test.exs @@ -278,7 +278,16 @@ defmodule Scholar.Decomposition.KernelPCATest do x = Nx.as_type(x(), :f64) model = KernelPCA.fit(x, num_components: 2, kernel: :rbf) assert Nx.type(model.eigenvalues) == {:f, 64} - assert Nx.type(KernelPCA.transform(model, Nx.as_type(x_test(), :f64))) == {:f, 64} + + transformed = KernelPCA.transform(model, Nx.as_type(x_test(), :f64)) + assert Nx.type(transformed) == {:f, 64} + + # f64 gets close to the f64 SciPy/scikit-learn reference + assert_all_close( + transformed, + Nx.tensor([[-0.11224132, 0.01060868], [-0.49994934, -0.10212445]], type: :f64), + atol: 1.0e-5 + ) end test "kernel defaults to :linear" do @@ -291,8 +300,22 @@ defmodule Scholar.Decomposition.KernelPCATest do end test "num_components may equal num_samples" do - model = KernelPCA.fit(x(), num_components: 6, kernel: :linear) + # The centering forces one zero eigenvalue (up to floating-point noise, + # which can make it slightly negative), so this exercises the clipping + # in fit/2 and the zero-eigenvalue guard in transform/2. + model = KernelPCA.fit(x(), num_components: 6, kernel: :rbf) assert Nx.shape(model.eigenvalues) == {6} + assert Nx.to_number(model.eigenvalues[5]) == 0.0 + + z = KernelPCA.fit_transform(x(), num_components: 6, kernel: :rbf) + assert Nx.to_number(Nx.any(Nx.is_nan(z))) == 0 + + zt = KernelPCA.transform(model, x_test()) + assert Nx.to_number(Nx.any(Nx.is_nan(zt))) == 0 + assert Nx.to_number(Nx.any(Nx.is_infinity(zt))) == 0 + + # the zero-eigenvalue component projects to zero, as in scikit-learn + assert zt[[.., 5]] == Nx.tensor([0.0, 0.0]) end test "eigenvalues are sorted in decreasing order" do From d295550235312bde6d0d0db02defa143750bc3c7 Mon Sep 17 00:00:00 2001 From: Ricardo Carvalho Santos Date: Sat, 11 Jul 2026 21:58:13 -0300 Subject: [PATCH 4/5] Refine eigenvalues with the Rayleigh quotient Verifying against scikit-learn on a larger dataset (25 samples) exposed that Nx.LinAlg.eigh can stop before the eigenvalues converge: the leading eigenvalue was exact, but the following ones came out up to a few percent too small (their eigenpair residuals were around 1.0e-2 where the leading one was at 1.0e-6), throwing fit_transform/2 off by the same margin. The 6-sample data used by the other tests was too small to show this. The eigenvectors themselves were accurate, so the eigenvalues are now recomputed as the Rayleigh quotient of the renormalized eigenvectors, whose error is quadratic in the eigenvector error. On the 25-sample data this brings the eigenvalues from an absolute error of 2.7 (poly kernel) to 1.0e-10 against scikit-learn, and fit_transform from 0.3 to 1.0e-5. Adds a regression test with the 25-sample dataset and refreshes the doctest values, which moved by a few ULPs. --- lib/scholar/decomposition/kernel_pca.ex | 20 ++++-- .../scholar/decomposition/kernel_pca_test.exs | 66 +++++++++++++++++++ 2 files changed, 79 insertions(+), 7 deletions(-) diff --git a/lib/scholar/decomposition/kernel_pca.ex b/lib/scholar/decomposition/kernel_pca.ex index ab104b41..cbfebfc4 100644 --- a/lib/scholar/decomposition/kernel_pca.ex +++ b/lib/scholar/decomposition/kernel_pca.ex @@ -92,7 +92,7 @@ defmodule Scholar.Decomposition.KernelPCA do iex> x = Nx.tensor([[0.5, 0.2, 0.8], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) iex> kpca = Scholar.Decomposition.KernelPCA.fit(x, num_components: 2, kernel: :rbf) iex> kpca.eigenvalues - Nx.tensor([0.3644832372665405, 0.26074573397636414]) + Nx.tensor([0.3644832372665405, 0.26074567437171936]) """ deftransform fit(x, opts \\ []) do opts = NimbleOptions.validate!(opts, @opts_schema) @@ -161,7 +161,7 @@ defmodule Scholar.Decomposition.KernelPCA do iex> Scholar.Decomposition.KernelPCA.transform(kpca, Nx.tensor([[0.5, 0.5, 0.5]])) Nx.tensor( [ - [0.12500189244747162, 0.029097510501742363] + [0.12500189244747162, 0.029097504913806915] ] ) """ @@ -213,10 +213,10 @@ defmodule Scholar.Decomposition.KernelPCA do iex> Scholar.Decomposition.KernelPCA.fit_transform(x, num_components: 2, kernel: :rbf) Nx.tensor( [ - [-0.13561560213565826, -0.16519638895988464], - [0.021600963547825813, 0.44114428758621216], - [0.4687873125076294, -0.15764620900154114], - [-0.35477250814437866, -0.11830171197652817] + [-0.13561560213565826, -0.16519635915756226], + [0.021600963547825813, 0.4411441683769226], + [0.4687873125076294, -0.15764616429805756], + [-0.35477250814437866, -0.11830168217420578] ] ) """ @@ -283,9 +283,15 @@ defmodule Scholar.Decomposition.KernelPCA do {eigenvalues, eigenvectors} = Nx.LinAlg.eigh(symmetric_kernel, eps: eigh_eps) order = Nx.argsort(eigenvalues, direction: :desc) - eigenvalues = Nx.take(eigenvalues, order)[0..(num_components - 1)] eigenvectors = Nx.take(eigenvectors, order, axis: 1)[[.., 0..(num_components - 1)]] + # eigh may stop before the eigenvalues fully converge, while the + # eigenvectors are already accurate, so the eigenvalues are recomputed as + # the Rayleigh quotient of the (renormalized) eigenvectors, whose error + # is quadratic in the eigenvector error + eigenvectors = eigenvectors / Nx.LinAlg.norm(eigenvectors, axes: [0]) + eigenvalues = Nx.sum(eigenvectors * Nx.dot(symmetric_kernel, eigenvectors), axes: [0]) + # negative eigenvalues carry no usable projection, either floating-point # noise around zero (the centering itself forces one zero eigenvalue) or # produced by an indefinite kernel such as sigmoid, so they are clipped diff --git a/test/scholar/decomposition/kernel_pca_test.exs b/test/scholar/decomposition/kernel_pca_test.exs index 15f1e9e7..bd14910c 100644 --- a/test/scholar/decomposition/kernel_pca_test.exs +++ b/test/scholar/decomposition/kernel_pca_test.exs @@ -323,6 +323,72 @@ defmodule Scholar.Decomposition.KernelPCATest do sorted = model.eigenvalues |> Nx.to_flat_list() |> Enum.sort(:desc) assert Nx.to_flat_list(model.eigenvalues) == sorted end + + test "matches scikit-learn on a larger dataset" do + # On more samples eigh stops before the eigenvalues converge, which the + # Rayleigh quotient refinement compensates for; without it the second + # eigenvalue here comes out around 20.73. + x = + Nx.tensor( + [ + [0.1236, 1.8521, 1.196, 0.796], + [-0.5319, -0.532, -0.8257, 1.5985], + [0.8033, 1.1242, -0.9382, 1.9097], + [1.4973, -0.363, -0.4545, -0.4498], + [-0.0873, 0.5743, 0.2958, -0.1263], + [0.8356, -0.5815, -0.1236, 0.0991], + [0.3682, 1.3555, -0.401, 0.5427], + [0.7772, -0.8606, 0.8226, -0.4884], + [-0.8048, 1.8467, 1.8969, 1.4252], + [-0.0862, -0.707, 1.0527, 0.3205], + [-0.6339, 0.4855, -0.8968, 1.728], + [-0.2237, 0.9876, -0.0649, 0.5602], + [0.6401, -0.4454, 1.9088, 1.3254], + [1.8185, 1.6845, 0.7937, 1.7656], + [-0.7345, -0.4121, -0.8643, -0.024], + [0.166, -0.186, 1.4862, 0.0703], + [-0.1572, 0.6281, -0.5772, 1.4066], + [-0.7763, 1.9607, 1.3167, -0.4039], + [-0.9834, 1.4464, 1.1206, 1.187], + [1.3138, -0.7779, 0.0754, -0.6524], + [1.5893, 0.8699, -0.0073, -0.8093], + [-0.0671, -0.0245, 1.1888, 0.9127], + [1.6616, 0.4166, -0.6412, 1.1397], + [1.2824, 0.6838, 1.3129, 0.4814], + [0.5682, 0.2826, -0.9237, -0.6763] + ], + type: :f64 + ) + + x_test = + Nx.tensor( + [ + [-0.9057, 0.9092, -0.0569, 0.5257], + [1.7227, -0.2521, 0.2311, 1.2667], + [-0.3136, -0.7691, -0.1307, -0.5163] + ], + type: :f64 + ) + + model = KernelPCA.fit(x, num_components: 4, kernel: :linear) + + assert_all_close( + model.eigenvalues, + Nx.tensor([30.12812417, 20.80703667, 15.30404198, 11.71610252], type: :f64) + ) + + assert_all_close( + KernelPCA.transform(model, x_test), + Nx.tensor( + [ + [0.62235235, 0.56524336, -0.82260582, -0.7056329], + [-0.72708075, 0.09505719, 0.88921225, 1.27344856], + [-1.18468312, -0.25757465, -1.27491476, -0.37172622] + ], + type: :f64 + ) + ) + end end describe "input validation" do From fd1ffac556fb3ff56b550d1e871f18a0294d293b Mon Sep 17 00:00:00 2001 From: Ricardo Carvalho Santos Date: Sat, 11 Jul 2026 23:11:53 -0300 Subject: [PATCH 5/5] Keep eigenvalues sorted after the Rayleigh refinement The component order comes from sorting the raw eigh eigenvalues, but the final values are the refined ones, so two eigenvalues that eigh resolved within its error could end up out of decreasing order, breaking the documented invariant. They are sorted again after the refinement. Also documents the kernel option fields stored in the struct. --- lib/scholar/decomposition/kernel_pca.ex | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/lib/scholar/decomposition/kernel_pca.ex b/lib/scholar/decomposition/kernel_pca.ex index cbfebfc4..fbd16bfe 100644 --- a/lib/scholar/decomposition/kernel_pca.ex +++ b/lib/scholar/decomposition/kernel_pca.ex @@ -87,6 +87,9 @@ defmodule Scholar.Decomposition.KernelPCA do * `:kernel_fit_all` - Mean of the whole training kernel matrix. + * `:kernel`, `:gamma`, `:degree`, `:coef0` - The kernel options used to fit + the model, applied again when computing the kernel of new samples. + ## Examples iex> x = Nx.tensor([[0.5, 0.2, 0.8], [1.0, 0.5, 0.2], [0.3, 1.0, 0.7], [0.9, 0.1, 1.0]]) @@ -292,6 +295,12 @@ defmodule Scholar.Decomposition.KernelPCA do eigenvectors = eigenvectors / Nx.LinAlg.norm(eigenvectors, axes: [0]) eigenvalues = Nx.sum(eigenvectors * Nx.dot(symmetric_kernel, eigenvectors), axes: [0]) + # the refinement can reorder eigenvalues that eigh resolved within its + # error, so they are sorted again to keep the decreasing order + order = Nx.argsort(eigenvalues, direction: :desc) + eigenvalues = Nx.take(eigenvalues, order) + eigenvectors = Nx.take(eigenvectors, order, axis: 1) + # negative eigenvalues carry no usable projection, either floating-point # noise around zero (the centering itself forces one zero eigenvalue) or # produced by an indefinite kernel such as sigmoid, so they are clipped