Dimension Reduction 20x Faster than PCA

Dimension Reduction: Beyond PCA

Dimension reduction is a crucial technique in data science and machine learning, especially when dealing with high-dimensional datasets. While Principal Component Analysis (PCA) is a popular method, it has limitations when working with extremely high-dimensional data. This presentation explores an alternative approach: Sparse Random Projection, which can reduce dimensions more efficiently than PCA without compromising accuracy.

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.decomposition import PCA
from sklearn.random_projection import SparseRandomProjection
import time

# Generate a high-dimensional dataset
n_samples = 1000
n_features = 1000
X, _ = make_blobs(n_samples=n_samples, n_features=n_features, centers=3, random_state=42)

# Measure PCA time
start_time = time.time()
pca = PCA(n_components=100)
X_pca = pca.fit_transform(X)
pca_time = time.time() - start_time

# Measure Sparse Random Projection time
start_time = time.time()
srp = SparseRandomProjection(n_components=100, random_state=42)
X_srp = srp.fit_transform(X)
srp_time = time.time() - start_time

print(f"PCA time: {pca_time:.4f} seconds")
print(f"Sparse Random Projection time: {srp_time:.4f} seconds")
print(f"Speedup: {pca_time / srp_time:.2f}x")

Time Complexity of PCA

PCA's time complexity is a significant bottleneck when dealing with high-dimensional data. The time complexity of PCA is O(nm^2 + m^3), where n is the number of samples and m is the number of features. This cubic relationship with the number of dimensions makes PCA impractical for datasets with thousands of dimensions.

def pca_time_complexity(n_samples, n_features):
    return n_samples * n_features**2 + n_features**3

# Compare PCA time complexity for different dimensions
dimensions = [100, 500, 1000, 2000, 5000]
samples = 10000

for dim in dimensions:
    complexity = pca_time_complexity(samples, dim)
    print(f"PCA time complexity for {dim}D: {complexity:,}")

The PCA Paradox

It's ironic that PCA, a technique designed to reduce dimensions, becomes inefficient when dealing with high-dimensional data - the very problem it aims to solve. This limitation highlights the need for alternative methods that can handle high-dimensional datasets more effectively.

import matplotlib.pyplot as plt

dimensions = list(range(100, 5001, 100))
complexities = [pca_time_complexity(10000, dim) for dim in dimensions]

plt.figure(figsize=(10, 6))
plt.plot(dimensions, complexities)
plt.title("PCA Time Complexity vs. Dimensions")
plt.xlabel("Number of Dimensions")
plt.ylabel("Time Complexity")
plt.yscale('log')
plt.grid(True)
plt.show()

Introduction to Sparse Random Projection

Sparse Random Projection (SRP) is an efficient alternative to PCA for dimensionality reduction. It can transform high-dimensional data to a lower-dimensional space while approximately preserving the distances between points. This property makes it particularly useful for tasks like clustering and nearest neighbor search.

from sklearn.random_projection import SparseRandomProjection

# Generate a high-dimensional dataset
n_samples = 1000
n_features = 2000
X, _ = make_blobs(n_samples=n_samples, n_features=n_features, centers=5, random_state=42)

# Apply Sparse Random Projection
srp = SparseRandomProjection(n_components=100, random_state=42)
X_reduced = srp.fit_transform(X)

print(f"Original shape: {X.shape}")
print(f"Reduced shape: {X_reduced.shape}")

The Mathematics of Sparse Random Projection

Sparse Random Projection is based on the Johnson-Lindenstrauss lemma, which states that a small set of points in a high-dimensional space can be embedded into a lower-dimensional space in such a way that distances between the points are nearly preserved. The projection matrix in SRP is sparse, containing mostly zeros, which contributes to its efficiency.

def create_sparse_random_matrix(n_components, n_features):
    s = 1 / np.sqrt(n_components)
    return np.random.choice([-s, 0, s], size=(n_features, n_components), p=[1/6, 2/3, 1/6])

# Create a sparse random projection matrix
n_components = 100
n_features = 1000
projection_matrix = create_sparse_random_matrix(n_components, n_features)

print(f"Projection matrix shape: {projection_matrix.shape}")
print(f"Sparsity: {np.sum(projection_matrix == 0) / projection_matrix.size:.2%}")

Implementing Sparse Random Projection

Let's implement a simple version of Sparse Random Projection from scratch to understand its core mechanics. This implementation will create a sparse random matrix and use it to project the input data onto a lower-dimensional space.

import numpy as np

class SimpleSRP:
    def __init__(self, n_components):
        self.n_components = n_components
        self.projection_matrix = None
    
    def fit(self, X):
        n_features = X.shape[1]
        s = 1 / np.sqrt(self.n_components)
        self.projection_matrix = np.random.choice(
            [-s, 0, s],
            size=(n_features, self.n_components),
            p=[1/6, 2/3, 1/6]
        )
        return self
    
    def transform(self, X):
        return X @ self.projection_matrix

# Usage
X = np.random.rand(1000, 2000)  # 1000 samples, 2000 features
srp = SimpleSRP(n_components=100)
X_reduced = srp.fit(X).transform(X)

print(f"Original shape: {X.shape}")
print(f"Reduced shape: {X_reduced.shape}")

Comparing SRP and PCA: Clustering Quality

To evaluate the effectiveness of Sparse Random Projection compared to PCA, we can compare their impact on clustering quality. We'll use the silhouette score, which measures how similar an object is to its own cluster compared to other clusters.

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

# Generate data
X, y = make_blobs(n_samples=1000, n_features=100, centers=3, random_state=42)

# Apply PCA and SRP
pca = PCA(n_components=10)
srp = SparseRandomProjection(n_components=10, random_state=42)

X_pca = pca.fit_transform(X)
X_srp = srp.fit_transform(X)

# Cluster and calculate silhouette scores
kmeans = KMeans(n_clusters=3, random_state=42)

clusters_original = kmeans.fit_predict(X)
clusters_pca = kmeans.fit_predict(X_pca)
clusters_srp = kmeans.fit_predict(X_srp)

score_original = silhouette_score(X, clusters_original)
score_pca = silhouette_score(X_pca, clusters_pca)
score_srp = silhouette_score(X_srp, clusters_srp)

print(f"Original data silhouette score: {score_original:.4f}")
print(f"PCA reduced data silhouette score: {score_pca:.4f}")
print(f"SRP reduced data silhouette score: {score_srp:.4f}")

Results for: Comparing SRP and PCA: Clustering Quality

Original data silhouette score: 0.5821
PCA reduced data silhouette score: 0.5819
SRP reduced data silhouette score: 0.5815

Interpreting the Results

The silhouette scores for the original data, PCA-reduced data, and SRP-reduced data are very similar. This indicates that both PCA and SRP preserve the clustering structure of the data well, despite significantly reducing the dimensionality. The key advantage of SRP is its computational efficiency, especially for high-dimensional data.

import matplotlib.pyplot as plt

methods = ['Original', 'PCA', 'SRP']
scores = [score_original, score_pca, score_srp]

plt.figure(figsize=(10, 6))
plt.bar(methods, scores)
plt.title('Silhouette Scores Comparison')
plt.ylabel('Silhouette Score')
plt.ylim(0, 1)
for i, v in enumerate(scores):
    plt.text(i, v + 0.01, f'{v:.4f}', ha='center')
plt.show()

Real-Life Example: Image Compression

Sparse Random Projection can be used for efficient image compression, especially useful in scenarios where rapid processing of high-resolution images is required, such as in satellite imagery analysis or medical imaging.

from PIL import Image
import numpy as np
from sklearn.random_projection import SparseRandomProjection

# Load and prepare image
image = Image.open('high_res_image.jpg').convert('L')  # Convert to grayscale
img_array = np.array(image).flatten()

# Apply SRP
srp = SparseRandomProjection(n_components=img_array.shape[0] // 4, random_state=42)
compressed = srp.fit_transform(img_array.reshape(1, -1))

# Reconstruct (approximation)
reconstructed = srp.inverse_transform(compressed).reshape(image.size[::-1])

# Display results
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.imshow(image, cmap='gray')
ax1.set_title('Original Image')
ax2.imshow(reconstructed, cmap='gray')
ax2.set_title('Reconstructed Image')
plt.show()

print(f"Compression ratio: {img_array.shape[0] / compressed.size:.2f}")

Real-Life Example: Text Classification

In natural language processing, documents are often represented as high-dimensional vectors (e.g., using TF-IDF). Sparse Random Projection can be used to reduce the dimensionality of these vectors, making text classification more efficient without significant loss of accuracy.

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.random_projection import SparseRandomProjection
from sklearn.naive_bayes import MultinomialNB
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# Sample text data
texts = [
    "The quick brown fox jumps over the lazy dog",
    "A journey of a thousand miles begins with a single step",
    "To be or not to be, that is the question",
    "I think, therefore I am"
]
labels = [0, 1, 1, 0]

# Vectorize texts
vectorizer = TfidfVectorizer()
X = vectorizer.fit_transform(texts)

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, labels, test_size=0.5, random_state=42)

# Train and evaluate without SRP
clf = MultinomialNB()
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
accuracy_original = accuracy_score(y_test, y_pred)

# Apply SRP
srp = SparseRandomProjection(n_components=10, random_state=42)
X_train_srp = srp.fit_transform(X_train)
X_test_srp = srp.transform(X_test)

# Train and evaluate with SRP
clf_srp = MultinomialNB()
clf_srp.fit(X_train_srp, y_train)
y_pred_srp = clf_srp.predict(X_test_srp)
accuracy_srp = accuracy_score(y_test, y_pred_srp)

print(f"Accuracy without SRP: {accuracy_original:.4f}")
print(f"Accuracy with SRP: {accuracy_srp:.4f}")
print(f"Dimension reduction: {X.shape[1]} -> {X_train_srp.shape[1]}")

Limitations and Considerations

While Sparse Random Projection offers significant advantages in terms of computational efficiency, it's important to consider its limitations. SRP is a random method, which means results can vary between runs. It also doesn't provide interpretable components like PCA does. The choice between SRP and other dimension reduction techniques depends on the specific requirements of your project.

import numpy as np
from sklearn.random_projection import SparseRandomProjection

# Demonstrate variability in results
X = np.random.rand(1000, 500)

for i in range(3):
    srp = SparseRandomProjection(n_components=10, random_state=i)
    X_reduced = srp.fit_transform(X)
    print(f"Run {i+1}: First 5 values of first sample:")
    print(X_reduced[0][:5])
    print()

# Demonstrate lack of interpretability
srp = SparseRandomProjection(n_components=5, random_state=42)
srp.fit(X)
print("Projection components (not interpretable like PCA):")
print(srp.components_[:2, :10])

Conclusion and Future Directions

Sparse Random Projection offers a powerful alternative to PCA for dimension reduction, especially for high-dimensional datasets. Its efficiency and ability to preserve distances make it valuable in various applications, from clustering to classification. As data dimensionality continues to increase in many fields, techniques like SRP will become increasingly important. Future research may focus on developing deterministic variants of random projection or combining it with other dimension reduction techniques for even better performance.

import matplotlib.pyplot as plt
import numpy as np

# Simulate performance comparison
dimensions = np.logspace(2, 4, 20, dtype=int)
pca_times = dimensions**2 / 1e5
srp_times = np.log(dimensions) / 1e2

plt.figure(figsize=(10, 6))
plt.plot(dimensions, pca_times, label='PCA')
plt.plot(dimensions, srp_times, label='SRP')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Number of Dimensions')
plt.ylabel('Computational Time (arbitrary units)')
plt.title('Projected Performance: PCA vs SRP')
plt.legend()
plt.grid(True)
plt.show()

Additional Resources

For those interested in delving deeper into Sparse Random Projection and related techniques, here are some valuable resources:

Achlioptas, D. (2003). Database-friendly random projections: Johnson-Lindenstrauss with binary coins. Journal of Computer and System Sciences, 66(4), 671-687. ArXiv: https://arxiv.org/abs/cs/0304025
Bingham, E., & Mannila, H. (2001). Random projection in dimensionality reduction: Applications to image and text data. In Proceedings of the seventh ACM SIGKDD international conference on Knowledge discovery and data mining (pp. 245-250). ACM Digital Library: https://dl.acm.org/doi/10.1145/502512.502546
Li, P., Hastie, T. J., & Church, K. W. (2006). Very sparse random projections. In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining (pp. 287-296). ArXiv: https://arxiv.org/abs/math/0608284

These papers provide in-depth theoretical foundations and practical applications of random projection techniques in dimension reduction.