I’m performing a bootstrapped principal factor analysis (PFA) on a dataset and encountering significant issues with aligning factors across bootstrap samples due to sign flipping and factor swapping, particularly with factors 3 and 4.

The Problem:

Sign Flipping: The signs of factor loadings arbitrarily flip between bootstrap samples, making it difficult to aggregate results.
Factor Swapping: Factors 3 and 4 often swap places between bootstraps, leading to mismatched assignments.
Aggregation Issues: These inconsistencies result in bimodal distributions of loading values when plotting histograms across bootstraps, indicating misalignment.
Approaches I’ve Tried:

Aligning Signs to a Reference:

Used the first bootstrap sample as a reference.
Flipped signs of factors in subsequent bootstraps if more than half of the signs differed.
Issue: Did not fully resolve the problem; factors were still misassigned.

Hungarian Algorithm with Correlation-Based Costs:

Computed absolute correlations between bootstrap factors and reference factors.
Used the negative absolute correlations as costs for the assignment problem.
Adjusted signs based on the sign of the correlation.
Issue: Factors 3 and 4 still got mixed up across bootstraps.

Clustering Factors:

Collected all factors from all bootstraps.
Standardized and clustered them using K-Means clustering.
Reordered factors based on cluster assignments.
Adjusted signs to align with cluster centers.
Issue: Factors remained misaligned; clustering did not produce consistent factor assignments.
Code Snippet:

Here’s the latest version of my post_process function using the clustering approach:

def post_process(results, sample_id, data):
    import numpy as np
    import pandas as pd
    from sklearn.cluster import KMeans
    from sklearn.preprocessing import StandardScaler

    # Number of factors
    n_factors = results['loadings'][0].shape[1]

    # Collect all factors from all bootstraps
    all_loadings = []
    for idx in range(len(results['loadings'])):
        loadings = results['loadings'][idx].copy()

        # Ensure consistent sign (e.g., largest absolute loading is positive)
        for i in range(loadings.shape[1]):
            max_abs_index = np.argmax(np.abs(loadings[:, i]))
            if loadings[max_abs_index, i] < 0:
                loadings[:, i] *= -1

        # Transpose to have factors as rows
        all_loadings.append(loadings.T)

    # Stack all factors
    all_loadings = np.vstack(all_loadings)

    # Standardize the factors
    scaler = StandardScaler()
    all_loadings_scaled = scaler.fit_transform(all_loadings)

    # Perform clustering
    kmeans = KMeans(n_clusters=n_factors, random_state=0)
    cluster_labels = kmeans.fit_predict(all_loadings_scaled)

    # Assign cluster labels back to factors
    factor_idx = 0
    for idx in range(len(results['loadings'])):
        loadings = results['loadings'][idx]
        scores = results['scores'][idx]
        n_factors_in_sample = loadings.shape[1]

        # Get cluster labels for the current sample's factors
        cluster_ids = cluster_labels[factor_idx:factor_idx + n_factors_in_sample]
        factor_idx += n_factors_in_sample

        # Map clusters to consistent order
        unique_clusters = np.unique(cluster_ids)
        cluster_order = {cluster: i for i, cluster in enumerate(sorted(unique_clusters))}
        ordered_indices = np.array([cluster_order[cluster_id] for cluster_id in cluster_ids])

        # Reorder factors based on cluster labels
        order = np.argsort(ordered_indices)
        results['loadings'][idx] = loadings[:, order]
        results['scores'][idx] = scores[:, order]

        # Adjust signs to align with cluster centers
        for i in range(n_factors_in_sample):
            cluster_center = kmeans.cluster_centers_[cluster_ids[order[i]]]
            loading = results['loadings'][idx][:, i]
            if np.dot(loading, cluster_center) < 0:
                results['loadings'][idx][:, i] *= -1
                results['scores'][idx][:, i] *= -1

    # Aggregation code follows...
    # [Omitted for brevity]

What I’ve Observed:

Despite these efforts, factors 3 and 4 still appear to swap and misalign across bootstraps.
Histograms of loading values for these factors show bimodal distributions.
This suggests that the current methods are not effectively resolving the factor alignment and sign flipping issues.

My Question:

How can I effectively align factors across bootstrap samples in PFA to resolve sign flipping and factor swapping, especially when factors are similar and tend to mix (like factors 3 and 4 in my case)?

Are there alternative methods or best practices for handling factor alignment and sign flipping in bootstrapped factor analysis?
How can I ensure that each factor is consistently assigned and signed across bootstrap samples for accurate aggregation?
Is there a way to handle factors that are inherently unstable or highly correlated?
Any insights, suggestions, or examples of how to address this issue would be greatly appreciated!