Biplots in matrix format using pca

This is a snippet of my dataframe:

    species bill_length_mm  bill_depth_mm   flipper_length_mm     body_mass_g   predicted_species
0   Adelie       18                   18         181             3750                Chinstrap
1   Adelie       17                   17         186             3800                Adelie
2   Adelie       18                   18         195             3250                Gentoo
3   Adelie       0                    0           0               0                  Adelie
4   Chinstrap    19                   19         193             3450                Chinstrap
5   Chinstrap    20                   20         190             3650                Gentoo
6   Chinstrap    17                   17         181             3625                Adelie
7   Gentoo       19                   19         195             4675                Chinstrap
8   Gentoo       18                   18         193             3475                Gentoo
9   Gentoo       20                   20         190             4250                Gentoo

I want to make a biplot for my data, which would be something like this:

But I want to make a biplot for every species vs predicted_species matrix, so 9 subplots,same as above, I am not sure how that can be achieved. One way could be to split into dataframes, and make a biplot for each, but that isn't very efficient and difficult for comparison.

Can anyone provide some suggestions on how this could be done?

Combining the answer by Qiyun Zhu on how to plot a biplot with my answer on how to split the plot into the true vs. predicted subsets, you could do it like this:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# Load iris data.
iris = sns.load_dataset('iris')
X = iris.iloc[:, :4].values
y = iris.iloc[:, 4].values
features = iris.columns[:4]
targets = ['setosa', 'versicolor', 'virginica']

# Mock up some predictions.
iris['species_pred'] = (40 * ['setosa'] + 5 * ['versicolor'] + 5 * ['virginica']
                        + 40 * ['versicolor'] + 5 * ['setosa'] + 5 * ['virginica']
                        + 40 * ['virginica'] + 5 * ['versicolor'] + 5 * ['setosa'])

# Reduce features to two dimensions.
X_scaled = StandardScaler().fit_transform(X)
pca = PCA(n_components=2).fit(X_scaled)
X_reduced = pca.transform(X_scaled)
iris[['pc1', 'pc2']] = X_reduced


def biplot(x, y, data=None, **kwargs):
    # Plot data points.
    sns.scatterplot(data=data, x=x, y=y, **kwargs)
    
    # Calculate arrow parameters.
    loadings = pca.components_[:2].T
    pvars = pca.explained_variance_ratio_[:2] * 100
    arrows = loadings * np.ptp(X_reduced, axis=0)
    width = -0.0075 * np.min([np.subtract(*plt.xlim()), np.subtract(*plt.ylim())])

    # Plot arrows.
    horizontal_alignment = ['right', 'left', 'right', 'right']
    vertical_alignment = ['bottom', 'top', 'top', 'bottom']
    for (i, arrow), ha, va in zip(enumerate(arrows), 
                                  horizontal_alignment, vertical_alignment):
        plt.arrow(0, 0, *arrow, color='k', alpha=0.5, width=width, ec='none',
                  length_includes_head=True)
        plt.text(*(arrow * 1.05), features[i], ha=ha, va=va, 
                 fontsize='small', color='gray')

    
# Plot small multiples, corresponding to confusion matrix.
sns.set()
g = sns.FacetGrid(iris, row='species', col='species_pred', 
                  hue='species', margin_titles=True)
g.map(biplot, 'pc1', 'pc2')
plt.show()