Implement Principal Component Analysis (PCA) module

src/stats/CMakeLists.txt

@@ -1,5 +1,9 @@
 set(stats_cppFiles
     ../stdlib_linalg_constants.fypp
+    ../stdlib_sorting.fypp
+    ../stdlib_sorting_ord_sort.fypp
+    ../stdlib_sorting_sort_adjoint.fypp
+    ../stdlib_sorting_sort.fypp
     )
 
 set(stats_fppFiles
@@ -10,12 +14,15 @@ set(stats_fppFiles
     ../stdlib_error.fypp
     ../stdlib_linalg.fypp
     ../stdlib_linalg_diag.fypp
+    ../stdlib_linalg_eigenvalues.fypp
     ../stdlib_linalg_state.fypp
+    ../stdlib_linalg_svd.fypp
     ../stdlib_random.fypp
     ../stdlib_selection.fypp
     ../stdlib_string_type.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
@@ -29,4 +36,9 @@ set(stats_fppFiles
     stdlib_stats_var.fypp
 )
 
+set(f90Files
+    ../stdlib_sorting_radix_sort.f90
+    )
+
 configure_stdlib_target(stats "" stats_fppFiles stats_cppFiles)
+target_link_libraries(stats PUBLIC blas lapack $<$<BOOL:${STDLIB_BITSET}>:bitsets>)

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,67 @@ module stdlib_stats
     #:endfor
   end interface moment
 
+
+  #! PCA supports all real kinds available in stdlib's internal BLAS/LAPACK.
+  #! When WITH_XDP or WITH_QP are enabled, extended precision is also supported.
+  #:set PCA_KINDS_TYPES = [("sp", "real(sp)"), ("dp", "real(dp)")]
+  #:if WITH_XDP
+  #:set PCA_KINDS_TYPES = PCA_KINDS_TYPES + [("xdp", "real(xdp)")]
+  #:endif
+  #:if WITH_QP
+  #:set PCA_KINDS_TYPES = PCA_KINDS_TYPES + [("qp", "real(qp)")]
+  #:endif
+
+  interface pca
+    !! version: experimental
+    !!
+    !! Principal Component Analysis (PCA)
+    !! ([Specification](../page/specs/stdlib_stats.html#pca))
+    #:for k1, t1 in PCA_KINDS_TYPES
+      module subroutine pca_${k1}$(x, components, singular_values, x_mean, &
+                                  method, overwrite_x, err)
+        ${t1}$, intent(inout) :: x(:,:)
+        ${t1}$, intent(out) :: components(:,:)
+        real(${k1}$), 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 in PCA_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 in PCA_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,247 @@
+#:include "common.fypp"
+#! PCA supports all real kinds available in stdlib's internal BLAS/LAPACK.
+#! When WITH_XDP or WITH_QP are enabled, extended precision is also supported.
+#:set PCA_KINDS_TYPES = [("sp", "real(sp)"), ("dp", "real(dp)")]
+#:if WITH_XDP
+#:set PCA_KINDS_TYPES = PCA_KINDS_TYPES + [("xdp", "real(xdp)")]
+#:endif
+#:if WITH_QP
+#:set PCA_KINDS_TYPES = PCA_KINDS_TYPES + [("qp", "real(qp)")]
+#:endif
+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
+  use stdlib_sorting, only: sort_index
+  implicit none
+
+contains
+
+  ! Helper subroutine: Centers data in-place by subtracting the mean from each column
+  #:for k1, t1 in PCA_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 j = 1, n
+         do i = 1, m
+            x(i, j) = x(i, j) - mu(j)
+         end do
+      end do
+    end subroutine center_data_${k1}$
+  #:endfor
+
+  ! SVD-based PCA driver: computes principal components via SVD of centered data
+  #:for k1, t1 in PCA_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(:,:)
+      real(${k1}$), 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 in PCA_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(:,:)
+      real(${k1}$), intent(out) :: singular_values(:)
+      type(linalg_state_type), intent(out) :: err
+
+      integer(ilp) :: i, j, m
+      integer(ilp), allocatable :: idx(:)
+      ${t1}$ :: alpha, beta
+      real(${k1}$) :: scale_factor
+      real(${k1}$), allocatable :: lambda(:), lambda_copy(:)
+      ${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
+      scale_factor = 1.0_${k1}$ / real(max(n-1, 1), ${k1}$)
+      alpha = scale_factor
+      beta = zero
+      allocate(c(p, p))
+      c = zero
+      call syrk('U', 'T', p, n, alpha, x_centered, n, beta, c, p)
+
+      ! Fill lower triangle from upper triangle (syrk only fills upper)
+      do j = 1, p-1
+         do i = j+1, p
+            c(i, j) = c(j, i)
+         end do
+      end do
+
+      allocate(lambda(p))
+      allocate(vectors(p, p))
+      call eigh(c, lambda, vectors=vectors, err=err)
+
+      if (err%ok()) then
+         ! Sort eigenvalues in descending order
+         allocate(lambda_copy(p))
+         lambda_copy = -lambda
+         allocate(idx(p))
+         call sort_index(lambda_copy, idx)
+
+         ! Assign sorted results with safety bounds checks
+         m = min(size(components, 1, kind=ilp), p)
+         m = min(m, size(singular_values, 1, kind=ilp))
+         do i = 1, m
+            components(i, :) = vectors(:, idx(i))
+            if (lambda(idx(i)) > 0.0_${k1}$) then
+               singular_values(i) = sqrt(lambda(idx(i)) * real(n-1, ${k1}$))
+            else
+               singular_values(i) = 0.0_${k1}$
+            end if
+         end do
+      end if
+    end subroutine pca_eigh_driver_${k1}$
+  #:endfor
+
+  #:for k1, t1 in PCA_KINDS_TYPES
+    module subroutine pca_${k1}$(x, components, singular_values, x_mean, &
+                                method, overwrite_x, err)
+      ${t1}$, intent(inout) :: x(:,:)
+      ${t1}$, intent(out) :: components(:,:)
+      real(${k1}$), 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 in PCA_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 in PCA_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 j = 1, p
+            do i = 1, n
+               x_reconstructed(i, j) = x_reconstructed(i, j) + x_mean(j)
+            end do
+         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)