You are a professional machine learning specialist who can answer complex machine learning questions in a professional and clear way. I am doing a lab that is performing dimension reduction for MNIST dataset.I performed LDA on the dataset.

The code:
"import numpy as np
import matplotlib.pyplot as plt
class LDA:
    def __init__(self, n_components=None):
        self.n_components = n_components
        self.eig_vectors = None

    def fit_transform(self,X,y):        # to reduce the dimensions of the input data
        height, width = X.shape
        unique_classes = np.unique(y)
        num_classes = len(unique_classes)
        scatter_t = np.cov(X.T)*(height - 1)
        scatter_w = 0       #scatter_w matrix denotes the intra-class covariance

        for i in range(num_classes):
            class_items = np.flatnonzero(y == unique_classes[i])
            scatter_w = scatter_w + np.cov(X[class_items].T) * (len(class_items)-1)

        scatter_b = scatter_t - scatter_w       #scatter_b is the inter-class covariance matrix
        # An eigenvector is a vector whose direction remains unchanged when a linear transformation is applied to it.
        eig_values, eig_vectors = np.linalg.eigh(np.linalg.pinv(scatter_w).dot(scatter_b))
        pc = X.dot(eig_vectors[:,::-1][:,:self.n_components])       #The first n_components are selected using the slicing operation

# If n_components is equal to 2, we plot the two components, considering each vector as one axis
        if self.n_components == 2:
            if y is None:
                plt.scatter(pc[:,0],pc[:,1])
            else:
                colors = ['r','g','b']
                labels = np.unique(y)
                for color, label in zip(colors, labels):
                    class_data = pc[np.flatnonzero(y==label)]
                    plt.scatter(class_data[:,0],class_data[:,1],c=color)
            plt.show()
        return pc, eig_vectors, eig_values
# Load 100 training samples
from mnist import load_SampleMnist, load_mnist
#nSamples = 100
trainpath = '../../lab1/data/noisy_train_digits.mat'
testpath = '../../lab1/data/noisy_test_digits.mat'
train_data, train_labels, test_data, test_labels = load_mnist(trainpath,testpath) #,nSamples)
LDA_object = LDA(n_components=9)
Y, eig_vectors, eig_values = LDA_object.fit_transform(train_data,train_labels)
eig_vectors = eig_vectors[:, np.argsort(-eig_values)]
eig_vectors = eig_vectors[:, 0:9]
print(eig_vectors.shape)

Xhat = Y @ np.transpose(eig_vectors)
print(Xhat.shape)

# Plot the original and reconstructed digits
X = test_data
#np.random.shuffle(X)
fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, )
for i in range(5):
    img = X[i].reshape(28,28)
    ax[i].imshow(img, interpolation='nearest')
fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, )
for i in range(5):
    img = Xhat[i].reshape(28,28)
    ax[i].imshow(img, interpolation='nearest')

"

The projection and reconstruction processes are as follows:
w=x_pri=U^T (x-µ)
x_rec=Ux_prj+µ=Uw+µ, where μ is the global mean of facial images.
I want to investigate the effect of not performing mean subtraction on the projected images by setting the global mean to zero for reconstructed faces when setting µ=0 in the above equations. Show code.