A possibly incorrect implementation of RGS algorithm discussed in Owen, A. B. (1994). Controlling correlations in Latin hypercube samples. Journal of the American Statistical Association, 89(428), 1517-1522.

RGS.py

#!/usr/bin/env python
# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt


# centered case
def LHS(n, d):
    samples = np.tile(np.arange(n, dtype=np.float64), (d, 1)).reshape(d, n)
    for i in range(d):
        np.random.shuffle(samples[i, :])
    samples = samples.T
    samples += 0.5
    samples /= n
    return samples

def OALHS(n, d, n_iter=5):
    # Initialization with standard LHS
    samples = LHS(n, d)
    apply_RGS(samples, n_iter)
    return samples

def apply_RGS(samples, n_iter=3):
    samples = samples.T
    
    # Apply orthogonalization
    for _ in range(n_iter):
        forward(samples)
        backward(samples)
        
    samples = samples.T
    return samples

def forward(X):
    d, n = X.shape
    X[:, :] -= 0.5  # de-mean before orthogonalization
    for i in range(d):
        for j in range(i):
            # np.sum(X[j, :]**2) can be pre-computed as np.sum(((np.arange(n) + 0.5) / n - 0.5)**2)
            X[i, :] -= np.dot(X[i, :], X[j, :]) / np.sum(X[j, :]**2) * X[j, :]
        X[i, :] = (np.argsort(X[i, :]) + 0.5) / n - 0.5  # without this, orthogoanlization is perfect
    X[:, :] += 0.5  # restore the mean 
    #return X
    
def backward(X):
    d, n = X.shape
    X[:, :] -= 0.5
    for i in reversed(range(d)):
        for j in reversed(range(i+1, d)):
            X[i, :] -= np.dot(X[i, :], X[j, :]) / np.sum(X[j, :]**2) * X[j, :]
        X[i, :] = (np.argsort(X[i, :]) + 0.5) / n - 0.5  # without this, orthogoanlization is perfect
    X[:, :] += 0.5
    #return X

def plot(X, n, d):
    fig = plt.figure(figsize=[6,6])
    ax = fig.gca()
    ax.set_xticks(np.arange(0, 1, 1./n))
    ax.set_yticks(np.arange(0, 1., 1./n))

    plt.scatter(X[:, 0], X[:, 1])

    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.grid()


if __name__ == '__main__':
    n, d = 100, 5
    X_lhs = LHS(n, d)
    X_oalhs = apply_RGS(np.copy(X_lhs), n_iter=1)
    Corr_lhs = np.corrcoef(X_lhs.T)
    Corr_oalhs = np.corrcoef(X_oalhs.T)
    Frob_SS_oalhs = np.sum((Corr_oalhs - np.diag([1]*d))**2) / 2
    Frob_SS_lhs = np.sum((Corr_lhs - np.diag([1]*d))**2) / 2
    print("Correlation matrix for LHS:\n{}\n"
          "with sum of squares of lower triangular components: {}".format(
          Corr_lhs, Frob_SS_lhs))
    print()
    print("Correlation matrix for OA-LHS:\n{}\n"
          "with sum of squares of lower triangular components: {}".format(
          Corr_oalhs, Frob_SS_oalhs))

    # n, d = 10, 2
    # X = LHS(n, d)
    # plot(X, n, d)
    # plt.show()