Implement Principal Component Analysis (PCA) module

src/stats/CMakeLists.txt

@@ -5,6 +5,7 @@ set(stats_fppFiles
     stdlib_random.fypp
     stdlib_stats_corr.fypp
     stdlib_stats_cov.fypp
+    stdlib_stats_pca.fypp
     stdlib_stats_distribution_exponential.fypp
     stdlib_stats_distribution_normal.fypp
     stdlib_stats_distribution_uniform.fypp

src/stats/stdlib_stats.fypp

@@ -9,10 +9,12 @@ module stdlib_stats
   !! ([Specification](../page/specs/stdlib_stats.html))
   use stdlib_kinds, only: sp, dp, xdp, qp, &
       int8, int16, int32, int64
+  use stdlib_linalg_state, only: linalg_state_type
   implicit none
   private
   ! Public API
   public :: corr, cov, mean, median, moment, var
+  public :: pca, pca_transform, pca_inverse_transform
 
 
   interface corr
@@ -639,4 +641,57 @@ module stdlib_stats
     #:endfor
   end interface moment
 
+
+  interface pca
+    !! version: experimental
+    !!
+    !! Principal Component Analysis (PCA)
+    !! ([Specification](../page/specs/stdlib_stats.html#pca))
+    #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+      module subroutine pca_${k1}$(x, components, singular_values, x_mean, &
+                                  method, overwrite_x, err)
+        ${t1}$, intent(inout) :: x(:,:)
+        ${t1}$, intent(out) :: components(:,:)
+        ${t1}$, intent(out) :: singular_values(:)
+        ${t1}$, intent(out), optional :: x_mean(:)
+        character(*), intent(in), optional :: method
+        logical, intent(in), optional :: overwrite_x
+        type(linalg_state_type), intent(out), optional :: err
+      end subroutine pca_${k1}$
+    #:endfor
+  end interface pca
+
+
+  interface pca_transform
+    !! version: experimental
+    !!
+    !! Projects data into the reduced dimensional space
+    !! ([Specification](../page/specs/stdlib_stats.html#pca_transform))
+    #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+      module subroutine pca_transform_${k1}$(x, components, x_mean, x_transformed)
+        ${t1}$, intent(in) :: x(:,:)
+        ${t1}$, intent(in) :: components(:,:)
+        ${t1}$, intent(in), optional :: x_mean(:)
+        ${t1}$, intent(out) :: x_transformed(:,:)
+      end subroutine pca_transform_${k1}$
+    #:endfor
+  end interface pca_transform
+
+
+  interface pca_inverse_transform
+    !! version: experimental
+    !!
+    !! Reconstructs original data from the reduced space
+    !! ([Specification](../page/specs/stdlib_stats.html#pca_inverse_transform))
+    #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+      module subroutine pca_inverse_transform_${k1}$(x_reduced, components, x_mean, x_reconstructed)
+        ${t1}$, intent(in) :: x_reduced(:,:)
+        ${t1}$, intent(in) :: components(:,:)
+        ${t1}$, intent(in), optional :: x_mean(:)
+        ${t1}$, intent(out) :: x_reconstructed(:,:)
+      end subroutine pca_inverse_transform_${k1}$
+    #:endfor
+  end interface pca_inverse_transform
+
+
 end module stdlib_stats

src/stats/stdlib_stats_pca.fypp

@@ -0,0 +1,224 @@
+#:include "common.fypp"
+submodule (stdlib_stats) stdlib_stats_pca
+  use stdlib_kinds, only: sp, dp, xdp, qp
+  use stdlib_optval, only: optval
+  use stdlib_linalg, only: svd, eigh
+  use stdlib_linalg_constants, only: ilp
+  use stdlib_linalg_blas, only: gemm, syrk
+  use stdlib_linalg_state, only: linalg_state_type, LINALG_ERROR, LINALG_VALUE_ERROR
+  implicit none
+
+contains
+
+  ! Helper subroutine: Centers data in-place by subtracting the mean from each column
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    subroutine center_data_${k1}$(x, mu)
+      ${t1}$, intent(inout) :: x(:,:)
+      ${t1}$, intent(in) :: mu(:)
+      integer(ilp) :: i, j, m, n
+      m = size(x, 1, kind=ilp)
+      n = size(x, 2, kind=ilp)
+      do concurrent(i=1:m, j=1:n)
+         x(i, j) = x(i, j) - mu(j)
+      end do
+    end subroutine center_data_${k1}$
+  #:endfor
+
+  ! SVD-based PCA driver: computes principal components via SVD of centered data
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    subroutine pca_svd_driver_${k1}$(x_centered, n, p, components, singular_values, err)
+      use stdlib_blas_constants_${k1}$, only: zero
+      ${t1}$, intent(inout) :: x_centered(:,:)
+      integer(ilp), intent(in) :: n, p
+      ${t1}$, intent(out) :: components(:,:)
+      ${t1}$, intent(out) :: singular_values(:)
+      type(linalg_state_type), intent(out) :: err
+
+      integer(ilp) :: n_s, m
+      ${t1}$, allocatable :: s_tmp(:), vt_tmp(:,:)
+
+      ! Initialize outputs to zero to prevent uninitialized memory access
+      components = zero
+      singular_values = zero
+
+      n_s = min(n, p)
+      allocate(s_tmp(n_s))
+      allocate(vt_tmp(n_s, p))
+
+      call svd(x_centered, s_tmp, vt=vt_tmp, overwrite_a=.true., full_matrices=.false., err=err)
+
+      if (err%ok()) then
+         m = min(size(components, 1, kind=ilp), n_s)
+         components(1:m, :) = vt_tmp(1:m, :)
+         m = min(size(singular_values, 1, kind=ilp), n_s)
+         singular_values(1:m) = s_tmp(1:m)
+      end if
+    end subroutine pca_svd_driver_${k1}$
+  #:endfor
+
+  ! Eigendecomposition-based PCA driver: computes principal components via covariance matrix
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    subroutine pca_eigh_driver_${k1}$(x_centered, n, p, components, singular_values, err)
+      use stdlib_blas_constants_${k1}$, only: zero
+      ${t1}$, intent(in) :: x_centered(:,:)
+      integer(ilp), intent(in) :: n, p
+      ${t1}$, intent(out) :: components(:,:)
+      ${t1}$, intent(out) :: singular_values(:)
+      type(linalg_state_type), intent(out) :: err
+
+      integer(ilp) :: i, j, m
+      ${t1}$ :: alpha, beta
+      real(${k1}$), allocatable :: lambda(:)
+      ${t1}$, allocatable :: c(:,:), vectors(:,:)
+
+      ! Initialize outputs to zero to prevent uninitialized memory access
+      components = zero
+      singular_values = zero
+
+      ! Compute covariance matrix using BLAS syrk: C = (1/(n-1)) * X^T * X
+      alpha = 1.0_${k1}$ / real(max(n-1, 1), ${k1}$)
+      beta = zero
+      allocate(c(p, p), source=zero)
+      call syrk('U', 'T', p, n, alpha, x_centered, n, beta, c, p)
+
+      ! Fill lower triangle from upper triangle (syrk only fills upper)
+      allocate(lambda(p))
+      allocate(vectors(p, p))
+      call eigh(c, lambda, vectors=vectors, upper_a=.true., err=err)
+
+      if (err%ok()) then
+         ! LAPACK returns eigenvalues in ascending order.
+         ! Flip them to get descending order for PCA.
+         lambda = lambda(p:1:-1)
+         vectors = vectors(:, p:1:-1)
+
+         ! Assign results with safety bounds checks
+         m = min(size(components, 1, kind=ilp), p)
+         m = min(m, size(singular_values, 1, kind=ilp))
+
+         ! Components are eigenvectors as rows (transpose of vectors columns)
+         components(1:m, :) = transpose(vectors(:, 1:m))
+
+         ! Convert eigenvalues to singular values: s = sqrt(lambda * (n-1))
+         where (lambda(1:m) > 0.0_${k1}$)
+            singular_values(1:m) = sqrt(lambda(1:m) * real(n-1, ${k1}$))
+         elsewhere
+            singular_values(1:m) = 0.0_${k1}$
+         end where
+      end if
+    end subroutine pca_eigh_driver_${k1}$
+  #:endfor
+
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    module subroutine pca_${k1}$(x, components, singular_values, x_mean, &
+                                method, overwrite_x, err)
+      ${t1}$, intent(inout) :: x(:,:)
+      ${t1}$, intent(out) :: components(:,:)
+      ${t1}$, intent(out) :: singular_values(:)
+      ${t1}$, intent(out), optional :: x_mean(:)
+      character(*), intent(in), optional :: method
+      logical, intent(in), optional :: overwrite_x
+      type(linalg_state_type), intent(out), optional :: err
+
+      type(linalg_state_type) :: err0
+      integer(ilp) :: n, p
+      ${t1}$, allocatable :: mu(:), x_centered(:,:)
+      character(len=:), allocatable :: method_
+
+      n = size(x, 1, kind=ilp)
+      p = size(x, 2, kind=ilp)
+      method_ = trim(adjustl(optval(method, "svd")))
+
+      ! Calculate mean along dimension 1 (column means) using stdlib mean
+      allocate(mu(p))
+      mu = mean(x, 1)
+
+      ! Validate and assign x_mean if present
+      if (present(x_mean)) then
+         if (size(x_mean) < p) then
+            err0 = linalg_state_type("pca", LINALG_VALUE_ERROR, &
+                   "x_mean array has insufficient size:", size(x_mean), ", expected:", p)
+            call err0%handle(err)
+            return
+         end if
+         x_mean(1:p) = mu
+      end if
+
+      ! Method dispatch
+      select case (method_)
+      case ("svd")
+         if (optval(overwrite_x, .false.)) then
+            call center_data_${k1}$(x, mu)
+            call pca_svd_driver_${k1}$(x, n, p, components, singular_values, err0)
+         else
+            allocate(x_centered, source=x)
+            call center_data_${k1}$(x_centered, mu)
+            call pca_svd_driver_${k1}$(x_centered, n, p, components, singular_values, err0)
+         end if
+
+      case ("eig", "cov")
+         allocate(x_centered, source=x)
+         call center_data_${k1}$(x_centered, mu)
+         call pca_eigh_driver_${k1}$(x_centered, n, p, components, singular_values, err0)
+
+      case default
+         err0 = linalg_state_type("pca", LINALG_ERROR, "Unknown method: ", method_)
+      end select
+
+      ! Handle error state
+      call err0%handle(err)
+
+    end subroutine pca_${k1}$
+  #:endfor
+
+
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    module subroutine pca_transform_${k1}$(x, components, x_mean, x_transformed)
+      use stdlib_blas_constants_${k1}$, only: one, zero
+      ${t1}$, intent(in) :: x(:,:)
+      ${t1}$, intent(in) :: components(:,:)
+      ${t1}$, intent(in), optional :: x_mean(:)
+      ${t1}$, intent(out) :: x_transformed(:,:)
+
+      integer(ilp) :: n, p, nc
+      ${t1}$, allocatable :: x_centered(:,:)
+
+      n = size(x, 1, kind=ilp)
+      p = size(x, 2, kind=ilp)
+      nc = size(components, 1, kind=ilp)
+
+      allocate(x_centered, source=x)
+      if (present(x_mean)) call center_data_${k1}$(x_centered, x_mean)
+
+      ! x_transformed = x_centered * components^T using GEMM
+      call gemm('N', 'T', n, nc, p, one, x_centered, n, components, nc, zero, x_transformed, n)
+    end subroutine pca_transform_${k1}$
+  #:endfor
+
+
+  #:for k1, t1, ri, cpp in REAL_KINDS_TYPES
+    module subroutine pca_inverse_transform_${k1}$(x_reduced, components, x_mean, x_reconstructed)
+      use stdlib_blas_constants_${k1}$, only: one, zero
+      ${t1}$, intent(in) :: x_reduced(:,:)
+      ${t1}$, intent(in) :: components(:,:)
+      ${t1}$, intent(in), optional :: x_mean(:)
+      ${t1}$, intent(out) :: x_reconstructed(:,:)
+
+      integer(ilp) :: i, j, n, nc, p
+
+      n = size(x_reduced, 1, kind=ilp)
+      nc = size(x_reduced, 2, kind=ilp)
+      p = size(components, 2, kind=ilp)
+
+      ! x_reconstructed = x_reduced * components using GEMM
+      call gemm('N', 'N', n, p, nc, one, x_reduced, n, components, nc, zero, x_reconstructed, n)
+
+      if (present(x_mean)) then
+         do concurrent(i=1:n, j=1:p)
+            x_reconstructed(i, j) = x_reconstructed(i, j) + x_mean(j)
+         end do
+      end if
+    end subroutine pca_inverse_transform_${k1}$
+  #:endfor
+
+end submodule stdlib_stats_pca

test/stats/CMakeLists.txt

@@ -8,12 +8,14 @@ set(fppFiles
     test_distribution_uniform.fypp
     test_distribution_normal.fypp
     test_distribution_exponential.fypp
+    test_pca.fypp
 )
 
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
 ADDTEST(corr)
 ADDTEST(cov)
+ADDTEST(pca)
 ADDTEST(mean)
 ADDTEST(median)
 ADDTEST(moment)

test/stats/test_pca.fypp

@@ -0,0 +1,54 @@
+#:include "common.fypp"
+program test_pca
+    use stdlib_error, only: check
+    use stdlib_kinds, only: sp, dp, xdp, qp
+    use stdlib_stats, only: pca, pca_transform, pca_inverse_transform
+    use stdlib_linalg_state, only: linalg_state_type
+    implicit none
+
+    #:for k1 in REAL_KINDS
+    real(${k1}$), parameter :: ${k1}$tol = 1000 * epsilon(1._${k1}$)
+    #:endfor
+
+    #:for k1 in REAL_KINDS
+    call test_pca_${k1}$()
+    #:endfor
+
+contains
+
+    #:for k1 in REAL_KINDS
+    subroutine test_pca_${k1}$()
+        real(${k1}$) :: x(3, 2), components(2, 2), s(2), mu(2)
+        real(${k1}$) :: x_trans(3, 2), x_inv(3, 2)
+        type(linalg_state_type) :: err
+
+        ! Data: [1, 2], [3, 4], [5, 6]
+        x = reshape([1.0_${k1}$, 3.0_${k1}$, 5.0_${k1}$, 2.0_${k1}$, 4.0_${k1}$, 6.0_${k1}$], [3, 2])
+
+        ! Test SVD method
+        call pca(x, components, s, x_mean=mu, method="svd", err=err)
+        call check(err%ok(), "pca_${k1}$ svd err")
+        call check(all(abs(mu - [3.0_${k1}$, 4.0_${k1}$]) < ${k1}$tol), "pca_${k1}$ svd mean")
+        ! First component should be approx [0.707, 0.707] (or negative)
+        call check(abs(abs(components(1,1)) - 1.0_${k1}$/sqrt(2.0_${k1}$)) < ${k1}$tol, "pca_${k1}$ svd comp1")
+        call check(abs(s(1) - 4.0_${k1}$) < ${k1}$tol, "pca_${k1}$ svd s1")
+        call check(abs(s(2)) < ${k1}$tol, "pca_${k1}$ svd s2")
+
+        ! Test Transform
+        call pca_transform(x, components, mu, x_trans)
+        ! Second dimension should be zero
+        call check(all(abs(x_trans(:, 2)) < ${k1}$tol), "pca_${k1}$ transform")
+
+        ! Test Inverse Transform
+        call pca_inverse_transform(x_trans, components, mu, x_inv)
+        call check(all(abs(x_inv - x) < ${k1}$tol), "pca_${k1}$ inverse")
+
+        ! Test EIG method
+        call pca(x, components, s, method="eig", err=err)
+        call check(err%ok(), "pca_${k1}$ eig err")
+        call check(abs(s(1) - 4.0_${k1}$) < ${k1}$tol, "pca_${k1}$ eig s1")
+
+    end subroutine test_pca_${k1}$
+    #:endfor
+
+end program test_pca