diff --git a/Project.toml b/Project.toml
index 3ab2cb7..724a4b7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,20 +7,24 @@ version = "0.1.3"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 PEG = "12d937ae-5f68-53be-93c9-3a6f997a20a8"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Aqua = "0.6"
 Dictionaries = "0.3"
 DocStringExtensions = "0.9"
 LightXML = "0.9"
+LinearAlgebra = "1.11.0"
 MacroTools = "0.5"
 PEG = "1.0"
 RecipesBase = "1.0"
 Reexport = "1.0"
+Statistics = "1.11.1"
 julia = "1.6"
 
 [extras]
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 61780bb..6dfd8ca 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -5,6 +5,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 ## Unreleased
 
 - ![](https://img.shields.io/badge/BREAKING-red.svg) use a different type parameter for each field of membership functions. This allows membership functions to support mixed-type inputs like `GaussianMF(0, 0.1)`.
+- ![](https://img.shields.io/badge/new%20feature-green.svg) added Gustafson-Kessel clustering algorithm.
 
 ## v0.1.3 -- 2024-09-18
 
diff --git a/src/FuzzyLogic.jl b/src/FuzzyLogic.jl
index 1538858..5cfd177 100644
--- a/src/FuzzyLogic.jl
+++ b/src/FuzzyLogic.jl
@@ -1,6 +1,6 @@
 module FuzzyLogic
 
-using Dictionaries, Reexport
+using Dictionaries, Reexport, LinearAlgebra, Statistics
 
 include("docstrings.jl")
 include("intervals.jl")
@@ -28,7 +28,7 @@ export DifferenceSigmoidMF, LinearMF, GeneralizedBellMF, GaussianMF, ProductSigm
        KarnikMendelDefuzzifier, EKMDefuzzifier, IASCDefuzzifier, EIASCDefuzzifier,
        @mamfis, MamdaniFuzzySystem, @sugfis, SugenoFuzzySystem, set,
        LinearSugenoOutput, ConstantSugenoOutput,
-       fuzzy_cmeans,
+       fuzzy_cmeans, gustafson_kessel,
        compilefis,
        readfis
 
diff --git a/src/genfis.jl b/src/genfis.jl
index 30446de..4b8cf94 100644
--- a/src/genfis.jl
+++ b/src/genfis.jl
@@ -1,3 +1,4 @@
+
 """
 Performs fuzzy clustering on th data `X` using `N` clusters.
 
@@ -54,3 +55,102 @@ function fuzzy_cmeans(X::Matrix{T}, N::Int; m = 2.0, maxiter = 100,
     end
     return C, U
 end
+
+"""
+Performs fuzzy clustering on the data `X` using `C` clusters.
+
+### Input
+
+- `X` -- ``d × N`` matrix of data, each column is a data point
+- `C` -- number of clusters used.
+
+### Keyword arguments
+
+- `α` -- exponent of the fuzzy membership function, default `2.0`
+- `ρ` -- ``1 × C`` vector of normalization factors for cluster sizes, default `ones(C)`
+- `maxiter` -- maximum number of iterations, default `100`
+- `tol` -- absolute error for stopping condition. Stop if ``|Eₖ - Eₖ₊₁|≤tol``, where ``Eₖ``
+           is the cost function value at the ``k``:th iteration.
+- `stabilize` -- applies covariance matrix stabilization proposed by Babuska et. al., default `true`
+- `γ` -- regularization coefficient, only used when stabilize is true, default `0.1`
+- `β` -- condition number threshold, only used when stabilize is true, default `10^15`
+### Output
+
+- `V` -- ``d × N``matrix of centers, each column is the center of a cluster.
+- `W` -- ``N × C`` matrix of membership degrees, ``Wᵢⱼ`` tells has the membership degree of
+         the ``i``th point to the ``j``th cluster.
+"""
+function gustafson_kessel(X::Matrix{T}, C::Int; α=2.0, ρ = ones(C),
+    maxiter=100, tol=1e-5, stabilize=true, γ=0.1, β=10^15) where {T <: Real}
+
+    α > 1 || throw(ArgumentError("α must be greater than 1"))
+    size(ρ, 1) == C || throw(ArgumentError("length of ρ must equal to C"))
+    γ >= 0 && γ <= 1 || throw(ArgumentError("γ must be within the range [0, 1]"))
+
+    N = size(X, 2)
+    d = size(X, 1)
+    W = rand(float(T), N, C)
+    W ./= sum(W, dims=2)
+    distances = zeros(N, C)
+    e = 2 / (α - 1)
+    J = Inf
+    P0_det = stabilize ? det(cov(X')) ^(1/d) : zero(float(T))
+
+    V = Matrix{float(T)}(undef, d, C)
+
+    @inbounds for iter in 1:maxiter
+        pow_W = W .^ α
+        W_sum = sum(pow_W, dims=1)
+        V .= (X * pow_W) ./ W_sum
+        for (i, vi) in enumerate(eachcol(V))
+            P = zeros(d, d)
+            for (j, xj) in enumerate(eachcol(X))
+                sub = xj - vi
+                P += pow_W[j, i] .* sub * sub'
+            end
+
+            P ./= W_sum[1, i]
+            stabilize && _recalculate_cov_matrix(P, γ, β, P0_det)
+            A = (ρ[i] * det(P))^(1/d) * inv(P)
+
+            for (j, xj) in enumerate(eachcol(X))
+                sub = xj - vi
+                distances[j, i] = sub' * A * sub
+            end
+        end
+
+        J_new = 0
+        for (i, di) in enumerate(eachrow(distances))
+            idx = findfirst(==(0), di)
+
+            if idx === nothing
+                for (j, dij) in enumerate(di)
+                    denom = sum((dij ./ di[k]) .^ e for k in 1:C)
+                    W[i, j] = 1/denom
+                    J_new += W[i, j] * dij
+                end
+            else
+                result = zeros(N)
+                result[idx] = 1
+                W[i, :] .= result
+            end
+        end
+
+        abs(J - J_new) < tol && break
+        J = J_new
+    end
+
+    return V, W
+end
+
+function _recalculate_cov_matrix(P::Matrix{Float64}, γ::Number,
+                                β::Int, P0_det::Float64)
+    temp_P = (1 - γ) * P + γ * P0_det * I(size(P, 1))
+    eig = eigen(Symmetric(temp_P))
+    eig_vals = eig.values
+    max_val = maximum(eig_vals)
+
+    eig_vals[max_val ./ eig_vals .> β] .= max_val / β
+
+    P .= eig.vectors * diagm(eig_vals) * eig.vectors'
+end
diff --git a/test/test_genfis.jl b/test/test_genfis.jl
index 9838eea..f56e121 100644
--- a/test/test_genfis.jl
+++ b/test/test_genfis.jl
@@ -1,4 +1,8 @@
-using FuzzyLogic, Test
+using FuzzyLogic, Random, LinearAlgebra, Test
+
+function generate_cluster(μ::Vector{Float64}, σ::Matrix{Float64}, n::Int)
+    return cholesky(Symmetric(σ)).L * randn(2, n) .+ μ
+end
 
 @testset "fuzzy c-means" begin
     X = [-3 -3 -3 -2 -2 -2 -1 0 1 2 2 2 3 3 3;
@@ -7,3 +11,31 @@ using FuzzyLogic, Test
     @test sortslices(C; dims = 2)≈[-2.02767 2.02767; 0 0] atol=1e-3
     @test_throws ArgumentError fuzzy_cmeans(X, 3; m = 1)
 end
+
+@testset "Gustafson-Kessel" begin
+    Random.seed!(5)
+    μ1 = [-2; 1.0]
+    μ2 = [6.0; 3]
+    X1 = generate_cluster(μ1, diagm(fill(0.25, 2)), 200)
+
+    θ = π / 4
+    U = [cos(θ) -sin(θ); sin(θ) cos(θ)]
+    D = Diagonal([5, 0.25])
+
+    Σ = U * D * U'
+
+    X2 = generate_cluster(fill(0.0, 2), Σ, 300)
+    X3 = generate_cluster(-μ1, diagm(fill(0.5, 2)), 200)
+    X4 = generate_cluster(μ2, diagm(fill(0.1, 2)), 200)
+
+    X = hcat(X1, X2, X3, X4)
+
+    Random.seed!()
+    (C, U) = gustafson_kessel(X, 4; α = 2)
+
+    @test sortslices(
+        C, dims = 2)≈[-1.98861 -0.08146 2.00787 5.94698; 1.04687 -0.03456 -0.95882 3.02226] atol=1e-3
+    @test_throws ArgumentError gustafson_kessel(X, 4; α = 1)
+    @test_throws ArgumentError gustafson_kessel(X, 4; ρ = [1, 1, 1])
+    @test_throws ArgumentError gustafson_kessel(X, 4, γ = 2)
+end