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')

"

 When projecting facial images onto low-dimensional subspace, it is necessary to remove the global mean (mean face) of facial images, and during reconstruction, it necessary to add the global mean as an offset. Specifically, 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 display some reconstructed faces when setting µ=0 in the above equations (use the maximum number of eigenfaces). Show code.