Graph laplacian is familiar to computation science researcher, with which we can perform spectral analysis such as diffusion map, eigenmap or spectral clustering. Here, we discuss how to generalize the graph laplacian to it's high-order form, i.e., Hodge laplacian.
In the previous post we know that the graph laplacian can be obtianed by degree matrix \(D\) and adjacency matrix \(A\) such that: \[\begin{equation} L_0 = D - A \end{equation}\]
We can find the adjacency matrix and degree matrix of the graph blow
with nine vertices: 
such that:
\[\begin{equation} A = \begin{array}[r]{c |c c c c c c c c c } & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\\ \hline 1& 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2& 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 3& 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 4& 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0\\ 5& 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0\\ 6& 0 & 0 & 0 & 1 & 1 & 0 & 1 & 0 & 0\\ 7& 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 8& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 9& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ \end{array} \end{equation}\]
and
\[\begin{equation} D = \begin{array}[r]{c |c c c c c c c c c } & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\\ \hline 1& 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2& 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 3& 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 4& 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0\\ 5& 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0\\ 6& 0 & 0 & 0 & 0 & 0 & 3 & 0 & 0 & 0\\ 7& 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 8& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 9& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ \end{array} \end{equation}\]
And the graph laplacian thus can be calculated as:
\[\begin{equation} L_0 = D - A = \begin{array}[r]{c |c c c c c c c c c } & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\\ \hline 1& 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2& -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\ 3& 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 4& 0 & 0 & 0 & 2 & -1 & -1 & 0 & 0 & 0\\ 5& 0 & 0 & 0 & -1 & 2 & -1 & 0 & 0 & 0\\ 6& 0 & 0 & 0 & -1 & -1 & 3 & -1 & 0 & 0\\ 7& 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0\\ 8& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1\\ 9& 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1\\ \end{array} \end{equation}\]
However, there's anoth way to represent the graph laplacian. Before that, we first introduce the incidence matrix of a graph. First, we can convert an undirected graph to be directed just assign the direction of edges by book keeping order. Next, let's construct the incidence matrix, each row represent a vertex, and each colum is an edge. We set the entry to be 1 if the the edge enter the vertex -1 when leave the vertex. Set 0 if there's no connection between a vertex and an edge. Thus the incidence matrix of a graph \(G\) is \(B\) such that: \[\begin{equation} B_1 = \begin{array}[r]{c | c c c c c c c} & \cdots & [4,5] & [4,6]& [5,6]& [6,7]& \cdots\\ \hline \vdots & \cdots & & \cdots & & \cdots & \\ 4 & \cdots & -1 & -1 & 0 & 0 & \cdots\\ 5 & \cdots & 1 & 0 & -1 & 0 & \cdots\\ 6 & \cdots & 0 & 1 & 1 & -1 & \cdots\\ 7 & \cdots & 0 & 0 & 0 & 1 & \cdots\\ \vdots & \cdots & & \cdots & & \cdots & \\ \end{array} \end{equation}\] Next, let construct graph laplacian matrix using \(B_1\). Actually, the graph laplacian is： \[\begin{equation} L_0=B_1 B_1^\top = \begin{array}[r]{c |c c c c c c c c c } & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\\ \hline 1& 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2& -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\ 3& 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 4& 0 & 0 & 0 & 2 & -1 & -1 & 0 & 0 & 0\\ 5& 0 & 0 & 0 & -1 & 2 & -1 & 0 & 0 & 0\\ 6& 0 & 0 & 0 & -1 & -1 & 3 & -1 & 0 & 0\\ 7& 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0\\ 8& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1\\ 9& 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 1\\ \end{array} \end{equation}\] As we can see that, \(L_0\) is the same as the that obtained from degree matrix and adjacency matrix.
We know that for the application of diffusion map or spectral analysis, we have to drop the first eigenvector due to its same values and the eigenvaule is 0. The second eigenvector is also called fidler vector. However, these applications usually created a connected graph. Since the graph we show here is an unconnected graph. We can stop here to take a guess: Is the zero eigenvaule still there, or if so, how many zero eigenvaules?
Now, let's perform the eigenvaule decomposition that: \[\begin{equation} L_0 = Q\Lambda Q^\top \end{equation}\]
Interestingly, the top 3 eigenvalues are all zero and we also have 3
connected components in the graph. Is this a coincidence or a property
of graph laplacian decomposition? 
Furthermore, the top 3 eigenvectors corresponding to the zero
eigenvaules are also interesting. Notice that the first eigenvector can
differentiate the components \(\{4,5,6,7\}\) from the rest vectices.
Likewise, the second eigenvector select the connected component \(\{8,9\}\) and the third choose the
component \(\{1,2,3\}\). 
In fact, in the field of Topological Data Analysis (TDA), the \(L_0\) is the special case of hodge laplacian (\(L_k\)). The number of connected components is call 0-dimensional cycles. And the graph laplaican can capture these cycles. The 1-dimensional cycles are correspond to holes, we will go into detail in the next section.
We can see the graph laplacian zero-order of hodge laplacian, and the formalu can be represented as: \[\begin{equation} L_0 = \mathbf{0}^\top \mathbf{0} + B_1 B_1^\top \end{equation}\]
Similarly, we can obtian \(L_1\):
\[\begin{equation} L_1 = B_1^\top B_1 + B_2 B_2^\top \end{equation}\]
You must ask where is \(B_2\) coming from? We know that \(B_1\) captures the relationship between vertices and edges. Thus, \(B_2\) captures the relationship between edges and triangles. We can also define \(B_3\) to capture relationship between triangles and tetrahedron and so on and so forth. So what is a triangle or a tetrahedron in a graph? We would not go into the detail of the thory of Simplex and Simplicial complex. Here we just need to know that three connected vertices forms a triangle. Similarly, four connected vertices forms a tetrahedron which is a high-order of triangles. To define \(B_2\), of a graph \(G\), we would check the direction of an edge \(e_j\) to the triangle \(\bigtriangleup_q\) it beblongs, if it has the same direction as the triangle, the entry would be \(1\), if the direction is opposite, the entry would be -1, otherwise the entry would be zero. Specifically:
\[\begin{equation} {B}_2[j, q] = \begin{cases} 1 & \text{if } e_j \in \bigtriangleup_q \quad \text{with}\quad \text{same}\quad \text{direction} \\ -1 & \text{if } e_j \in \bigtriangleup_q \quad \text{with}\quad \text{opposite}\quad \text{direction} \\ 0 & \text{otherwise} \end{cases} \end{equation}\]
With the definition, we can obtian \(B_2\) of the graph aforementioned: \[\begin{equation} B_2 = \begin{array}[r]{c | c } & [4,5,6]\\ \hline \vdots & \cdots \\ [4,5] & 1\\ [4,6] & -1\\ [5,6]& 1\\ [6,7] & 0\\ \vdots & \cdots \\ \end{array} \end{equation}\]
We next introduce the normalized form and the decomposition of hodge 1-laplacian. The normalized form of hodge 1-laplacian is given:
\[\begin{equation} \mathcal{L}_1 = {D}_2 {B}_1^\top {D}_1^{-1} {B}_1 + {B}_2 {D}_3 {B}_2^\top {D}_2^{-1} \end{equation}\] where \(\mathbf{D}_1\) is the vertices degree matrix, \({D}_2\) is \(\max{(\text{diag}(|{B}_2| \mathbf{1}), \mathbf{I})}\) and \({D}_3\) is \(\frac{1}{3}\mathbf{I}\). Since the normalized \(L_1\) is not neccessarily symmetric, we next need to define the symmetric normalized Hodge 1-Laplacian such that: \[\begin{equation} \begin{aligned} \mathcal{L}_1^s & = {D}_2^{-1/2} \mathcal{L}_1 {D}_2^{1/2}\\ & = {D}_2^{1/2} {B}_1^\top {D}_1^{-1} {B}_1 {D}_2^{1/2} + {D}^{-1/2} {B}_2 {D}_3 {B}_2^\top {D}_2^{-1/2} \end{aligned} \end{equation}\]
We use the graph with three holes to present hodge 1-laplacian: 
We next can perform eigenvalues decomposition on \(\mathcal{L}_1\): \[\begin{equation}
\begin{aligned}
\mathcal{L}_1
& = \mathbf{D}_2^{1/2} \mathcal{L}_1^s
\mathbf{D}_2^{-1/2} \\
& = \mathbf{D}_2^{1/2} Q \Lambda Q^\top
\mathbf{D}_2^{-1/2} \\
& = \mathbf{U} \Lambda \mathbf{U}^{-1}
\end{aligned}
\end{equation}\] Interestingly, the top 3 eigenvector also all
zero which corresponding to the 3 holes, namely, the three 1-dimensional
cycles. 
When it comes to the eigenvectors, we can also notice that the top
three eigenvectors are around the three holes. 
In algebraic geometry, these the eigenvectors with zero eigenvaules are called harmonic function or harmonic. These harmonic function around holes is useful for some analysis like clustering to find the 1-dimensional cycles etc.
Today, we review the graph laplacian, and we can find the zeor eigenvalues and theirs corresponding eigenvectors can be used to find connected components. The high order graph laplacian hodge laplacian have the similar properties, we presented the hodge 1-laplaican and its eigenvalues decomposition. We can find that the zero eigenvalues indicates the number of holes of a graph. Furthermore, the corresponding eigenvectors with zero eigenvalues are around holes.
Eigenvalue decomposition of weighted graph laplacian is used to obtain diffusion map. The asymmetric matrix or transition matrix \(P=D^{-1}W\). The decomposition is: \[\begin{equation} P = D^{-1/2}S\Lambda S^\top D^{1/2} \end{equation}\] In fact, the random walk on the graph would give rise to the transition graph such that: \[\begin{equation} p^{t+1}= p^t D^{-1}W = p^t P \end{equation}\] where \(p\) is a initial state of the graph (vertices weight). That is, after a certain step of random walk on the graph, we would reach a steady state when any more random walk would not change the weight. Let \(Q=D^{-1/2}S\), then \(Q^{-1}=S^{\top} D^{1/2}\), we can derive \(P = Q\Lambda Q^{-1}\). The random walk above mentioned can then be represent: \[\begin{equation} \begin{aligned} P^t &= Q\Lambda Q^{-1}Q\Lambda Q^{-1}\cdots Q\Lambda Q^{-1} &= Q\Lambda^t Q \end{aligned} \end{equation}\] Since \(\Lambda\) is diagnoal, the random walk on a graph can be easily cacluted by using the eigenvalue decomposition outcome. The column of \(Q\) is the diffusion map dimensions.
To understand diffusion map, we here introduce an example data and run diffusion map on it. The data is a single cell fibroblasts to neuron data with 392 cells. We first download the data from NCBI and loaded the data to memory:
1 | import requests |
From the code as we can see that, we extract the cell types from the
data and stores it into the variable group. The rest part
of the data is the normalized gene expression count matrix.
To do the pre-processing, we first get the log-normalized count
matrix \(X\) and revert it to the
counts and apply log geomatric scaling to the count matrix as the
pre-processing then store it to matrix \(Y\). The geomatric scaling
implementation is: 1
2
3X = np.array(data.iloc[:, 5:]).T
X = np.power(2, X[np.apply_along_axis(np.var, 1, X)>0, :]) - 1
Y = np.log(gscale(X+0.5)).T After the pre-processing, we we next
run PCA to perform dimensionality reduction and use which to calculate
euclidean distances between cells, and then run diffusion map.1
2
3
4
5
6
7
8
9def gscale(X:np.ndarray) -> np.ndarray:
assert(X.all()>=0)
div_ = np.divide(X.T, np.apply_along_axis(lambda x:np.exp(np.mean(np.log(x))), 1,
X)).T
scale_ = np.apply_along_axis(np.median,0, div_)
sc = StandardScaler(with_mean=False)
sc.fit(X)
sc.scale_ = scale_
return sc.transform(X)
1 | pc = run_pca(Y, 100) |
We first take a look at the PCA figure which shows the first 2 PCs of
the dataset. PCA is a linear methods which cannot reflect cell fate
differentation process. 
Since the diffusion map applied the Gaussian kernal during the distances
calculation, it usually a better way to capture the cell differentation
events. Here we show diffusion dimension 1,2 and 1,3. It's clear that
1,2 captures the cell differentation from fibroblasts to neuron, and 1,3
captures the differentiation from firoblasts to myocytes. 
We can change the bandwidth from 7 to 100, which use the 100th
neareast neighbor as bandwidth instead of 7th. The following shows the
diffusion map that bandwidth=100. 
The eigenvalue decomposition of graph laplacian can also be used to
infer the pseudo time simulating the pseudo time of cell differentation
events. We here set top \(m\) cells
\(1/m\) of the progenitor cells (MEF)
and set other cell 0 as the initial state. Next, we perform random walk
to get new pesudotime \(u\) such that:
\[\begin{equation}
u = [\frac{1}{m}, \frac{1}{m}, \cdots, 0, \cdots 0](D^{-1}W)^t
\end{equation}\] By testing different number of random walk steps
we can check the new pseudo time \(u\).
Here we show time step equals to 1, 10, 100 and 200. From the figure we
can notice that after certain step, the pseudo time will not change
anymore. That means the random walk reaches the steady state. 
To test when can we reach the steady state, I use the graph we
mentioned in my last post Spectral
analysis (1):
Here we random generate two initial states (\(v_1\), \(v_2\)) to do the random walk such that: \[\begin{equation} \begin{aligned} v_1 &= [0.15,0.41,0.54,0.9,0.62,0.93,0.1,0.46,0.01,0.88]\\ v_2 &= [0.89,0.93,0.07,0.41,0.52,0.88,0.43,0.09,0.1,0.2] \end{aligned} \end{equation}\]
From the the figure below we can see that \(v_1\) reaches the steady state in 30 steps
whereas \(v_2\) reaches steady state in
100 steps. In total, all these two initial state will reach the steady
state that would not change any longer. 
The eigenvectors of graph Laplacian can also be used to do the
clustering. From the figure we can find clear cut using the 1st, 2nd and
the 3rd diffusion dimensions. In practice, we can use kmeans, DBSCAN,
leiden, louvain algorithm to perform clustering using the diffusion
dimensions. 
1 | def run_pca(mtx, n_components=2, random_state=2022): |
1 | def affinity(R, k=7, sigma=None, log=False): |
1 | def diffusionMaps(R,k=7,sigma=None): |
Newton used glass separted the white sunlight into red, orange,
yellow, green, blue, indigo and violet single colors. The light spectrum
analysis can help scientists to interpret the world. For example, we can
detect the elements of our solar system as well as far stars in the
universe. The spectrum analysis is also a field in mathmatjics. In the
graph thoery field, Laplacian matrix is used to represented a graph. We
can obtian features from the undirected graph below (wikipedia).
For example, we can check the degree of each vertex, which forms our Degree matrix \(D\) such that: \[\begin{equation} D = \begin{pmatrix} 2 & 0 & 0 & 0 & 0 & 0\\ 0 & 3 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 3 & 0 &0\\ 0 & 0 & 0 & 0 & 3 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \end{equation}\]
By checking the connection between all pairs nodes, we can create a Adjacency matrix:
\[\begin{equation} A = \begin{pmatrix} 0 & 1 & 0 & 0 & 1 & 0\\ 1 & 0 & 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 & 1 &1\\ 1 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 \end{pmatrix} \end{equation}\]
The graph Laplacian \(L\) extracts all useful information from a graph which is: \[\begin{equation} L=D-A =\begin{pmatrix} 2 & -1 & 0 & 0 & -1 & 0\\ -1 & 3 & -1 & 0 & -1 & 0\\ 0 & -1 & 2 & -1 & 0 & 0\\ 0 & 0 & -1 & 3 & -1 & -1\\ -1 & -1 & 0 & -1 & 3 & 0\\ 0 & 0 & 0 & -1 & 0 & 1 \end{pmatrix} \end{equation}\]
In fact, the graph Laplaican matrix is symmetric and also positive semidefinite (PSD), which means if we perform eigenvalue decomposition, the eigen values are all real and nonnegative. We can normlize the graph Laplaican by left multiplying the \(D^{-1}\) which will give rise to all \(1\)s of the diagnal entries of the matrix, namely: \[\begin{equation} \text{norm}({L}) = D^{-1}L = \begin{pmatrix}1 & -0.50 & 0 & 0 & -0.50 & 0\\ -0.33 & 1 & -0.33 & 0 & -0.33 & 0\\ 0 & -0.50 & 1 & -0.50 & 0 & 0\\ 0 & 0 & -0.33 & 1 & -0.33 &-0.33\\ -0.33 & -0.33 & 0 & -0.33 & 1 & 0\\ 0 & 0 & 0 & -1 & 0 & 1 \end{pmatrix} \end{equation}\] This matrix is called transition matrix. However, the matrix would not keep the symmetric and the PSD property after the normalization. We will come back to discuss the spectral for the normalized graph Laplacian.
In practice, graph are usually weighted. The weight between vertices
can be euclidean distance or other measures. The figure blow shows the
weights of the graph. Here we apply gussian kernel to the euclidean
distances between vertices such that: 
\[\begin{equation} w_{ij} =\exp (-r_{ij}^2/\sigma) = \exp \big(\frac{-\|x_i - x_j\|^2}{\sigma_i\sigma_j}\big) \end{equation}\] where \(r_{ij}\) is the euclidean distance between vetex \(i\) and \(j\), and \(sigma\) controls the bandwidth. We call the matrix \(W=(w_{ij})\) gaussian affinity matrix such that: \[\begin{equation} \begin{pmatrix} 0 & w_{12} & w_{13} & w_{14} & w_{15}& w_{16}\\ w_{21} & 0 & w_{23} & w_{24} & w_{25}& w_{26}\\ w_{31} & w_{32} & 0 & w_{34} & w_{35}& w_{36}\\ w_{41} & w_{42} & w_{43} & 0 & w_{45}& w_{46}\\ w_{51} & w_{52} & w_{53} & w_{54} & 0& w_{56}\\ w_{61} & w_{62} & w_{63} & w_{64} & w_{65}& 0 \end{pmatrix} \end{equation}\] The guassian kernel would enlarge the distance between too far vertices. Similar to unweighted matrix, we can also construct graph Laplaican matrix using the gaussian affinity matrix. First, we need find the weighted degree based on \(W\) such that: \[\begin{equation} d_{ii} = \sum_j{w_{ij}} \end{equation}\] With the diagnal degree matrix and affinity matrix, we now can have the weighted laplacian that: \[\begin{equation} L = D - W \end{equation}\] Likewise, we next give the normalized form of Laplacian such that: \[\begin{equation} \text{norm}{(L)}= D^{-1}L = I - D^{-1}W \end{equation}\] To facilitate the eigenvalue decomposition, we need apply trick to the asymmetric matrix \(D^{-1}L\). Since the eigenvectors of \(D^{-1}L\) and \(D^{-1}W\) are the same, we apply some trick to \(P = D^{-1}W\) to simplify the problem. Lets construct \(P'\) such that: \[\begin{equation} P' = D^{1/2} P D^{-1/2} = D^{-1/2}WD^{-1/2} \end{equation}\] It obvious \(P'\) is symmetric due to the symmetrisation of \(W\). We can perform eigenvalue decomposition on \(P'\) such that: \[\begin{equation} P' = S\Lambda S^\top \end{equation}\] Where S stores the eigenvectors of \(P'\) and the diagonals of \(\Lambda\) records the eigenvalues of \(P'\). We can also get the decompostion to \(P\) such that: \[\begin{equation} P = D^{-1/2}S\Lambda S^\top D^{1/2} \end{equation}\] Let \(Q=D^{-1/2}\), then \(Q^{-1}=S^{\top}D^{1/2}\). We therefore find the right and left eigenvector of \(P\) such that: \[\begin{equation} \psi = D^{-1/2} S \qquad \phi = S^{\top}D^{1/2} \end{equation}\] In fact, columns of \(\psi\) stores the spectral of the graph which also call diffusion map dimensions.
1
2
3
4
t = tensor([[1,2,4],[4,5,6]])
t.shape: (2,3)
t.ndim: 2
type: scalar, vector, matrix, tensor
1
2
3
4
5
random_tensor = torch.rand(size=(3, 4, 5))
zeros = torch.zeros(size=(3, 4))
ones = torch.ones(size=(3, 4))
zero_to_ten = torch.arange(start=0, end=10, step=1)
ten_zeros = torch.zeros_like(input=zero_to_ten) # same shape but all zeros
1
2
3
4
torch.float32/torch.float
torch.float16
torch.half
torch.float64/torch.double
1
2
device=’cuda’ if torch.cuda.is_available() else ‘cpu’
t = tensor([1,2,3], device=device)
1
2
tensor_A = torch.tensor([[1,2],[3,4],[5,6]],
dtype = torch.float32)
1
2
3
4
5
6
7
tensor = torch.tensor([1, 2, 3])
tensor + 10
torch.multiply(tensor, 10)
tensor * tensor # tensor([1,4,9])
tenorA @ tensorB # matrix multiplication -> tensor(14)
torch.matmul(tensor, tensor)/torch.mm # 1*1 + 2*2 + 3*3 = tensor(14)
tensor.T
1
2
3
4
5
6
7
8
#Torch.nn.Linear
y = x A^T + b
torch.manual_seed(42)
linear = torch.nn.Linear(in_features=2, # in_features = matches inner dimension of input
out_features=6) # out_features = describes outer value
x = tensor_A
output = linear(x)
x.shape, output, output.shape
1
2
3
4
5
6
7
8
9
10
11
12
tensor = torch.arange(10, 100, 10) # tensor([10, 20, 30, 40, 50, 60, 70, 80, 90])
tensor.argmax() # 8
tensor.argmin() # 0
tensor.type(torch.float16) # tensor([10., 20., 30., 40., 50., 60., 70., 80., 90.]
torch.reshape(new_shape) # -1 is to ask calculating automatically
tensor.view(new_shape) # return a new shape view
torch.stack(t, dim=0) # concate a sequence of tensors along a new dimension(dim)
torch.squeeze() # all into the first dimensions
torch.clamp() # min=min, max=max, limit the range
torch.unsqueeze()
torch.permute() # torch.Size([224, 224, 3]) -> torch.Size([3, 224, 224])
torch.permute_(), x.unsqueeze_() -> inplace operation
1
2
3
4
x[:, 0] # tensor([[1, 2, 3]])
x[:, :, 1] # tensor([[2,5,8]])
x[:, 1, 1] # tensor([5])
x[0, 0, :]/x[0][0] # tensor([1,2,3])
1
tensor = torch.from_numpy(array)
1
2
3
random seed
torch.manual_seed(seed=RANDOM_SEED)
torch.random.manual_seed(seed=RANDOM_SEED)
1
2
3
4
5
6
7
8
9
10
torch.autograd import Variable
.data, .grad, .grad_fn
x_tensor = torch.randn(10, 5)
y_tensor = torch.randn(10, 5)
x = Variable(x_tensor, requires_grad=True)
y = Variable(y_tensor, requires_grad=True)
z = torch.sum(x + y)
print(z.data) #-2.1379
print(z.grad_fn) #<SumBackward0 object at 0x10da636a0>
z.backward()
1
2
3
4
5
6
7
8
9
if torch.cuda.is_available():
device = "cuda" # Use NVIDIA GPU (if available)
elif torch.backends.mps.is_available():
device = "mps" # Use Apple Silicon GPU (if available)
else:
device = "cpu" # Default to CPU if no GPU is available
tensor.to(device)
tensor_on_gpu.cpu().numpy()
torch.nn
Contains all of the building blocks for computational graphs (essentially a series of computations executed in a particular way).
torch.nn.Parameter
Stores tensors that can be used with nn.Module. If requires_grad=True gradients (used for updating model parameters via gradient descent) are calculated automatically, this is often referred to as "autograd".
torch.nn.Module
The base class for all neural network modules, all the building blocks for neural networks are subclasses. If you're building a neural network in PyTorch, your models should subclass nn.Module. Requires a forward() method be implemented.
torch.optim
Contains various optimization algorithms (these tell the model parameters stored in nn.Parameter how to best change to improve gradient descent and in turn reduce the loss).
def forward()
All nn.Module subclasses require a forward() method, this defines the computation that will take place on the data passed to the particular nn.Module (e.g. the linear regression formula above).
1
2
3
4
5
6
7
8
9
10
11
12
class net(nn.Module):
__init__(self): super().__init__() ...
def forward(self, x: torch.Tensor) -> torch.Tensor:
return self.weights * x + self.bias
## expample
class LinearRegressModle(nn.Module):
def __init__(self):
super().__init__()
self.weights = nn.Parameter(torch.randn(1, required_grad=True, dtype=torch.float))
self.bias = nn.Parameter(torch.randn(1, required_grad=True, dtype=torch.float))
def forward(self, x:torch.Tensor) -> torch.Tensor:
return self.weight * x + self.bias
1
2
3
4
5
torch.manual_seed(42)
model_0 = LinearRegressionModel()
list(model_0.parameters()) # tensor([0.3367], requires_grad=True)
model_0.state_dict() # OrderedDict([('weights', tensor([0.3367])), ('bias', tensor([0.1288]))])
with torch.inference_mode(): y_preds = model_0(X_test) # run inference
1
2
loss_fn = nn.L1Loss() # MAE loss is same as L1Loss
optimizer = torch.optim.SGD(params=model_0.parameters(), lr=0.01) ## lr(learning rate)
| Step name | What does it do? | Code example | |
|---|---|---|---|
| 1 | Forward pass | The model goes through all of the training data once, performing its forward() function calculations. | model(x_train) |
| 2 | Calculate the loss | The model's outputs (predictions) are compared to the ground truth and evaluated to see how wrong they are. | loss = loss_fn(y_pred, y_train) |
| 3 | Zero gradients | The optimizers gradients are set to zero (they are accumulated by default) so they can be recalculated for the specific training step. | optimizer.zero_grad() |
| 4 | Perform backpropagation on the loss | Computes the gradient of the loss with respect for every model parameter to be updated (each parameter with requires_grad=True) | loss.backward() |
| 5 | Update the optimizer (gradient descent) | Update the parameters with requires_grad=True with respect to the loss gradients in order to improve them. | optimizer.step() |
1
2
3
4
5
6
7
for epoch in range(epoches):
model.train()
y_pred = model(X_train)
loss = loss_fn(y_pred, y_true)
optimizer.zero_grad()
loss.backward()
optimizer.step()
| Forward pass | The model goes through all of the training data once, performing its forward() function calculations. | model(x_test) | |
| Calculate the loss | The model's outputs (predictions) are compared to the ground truth and evaluated to see how wrong they are. | loss = loss_fn(y_pred, y_test) | |
| Calulate evaluation metrics (optional) | Alongisde the loss value you may want to calculate other evaluation metrics such as accuracy on the test set. | Custom functions |
Inferennce
model_0.eval() # Set the model in evaluation mode with torch.inference_mode(): y_preds = model_0(X_test)
torch.save
Saves a serialized object to disk using Python's pickle utility. Models, tensors and various other Python objects like dictionaries can be saved using torch.save.
torch.load
Uses pickle's unpickling features to deserialize and load pickled Python object files (like models, tensors or dictionaries) into memory. You can also set which device to load the object to (CPU, GPU etc).
torch.nn.Module.load_state_dict ## recommended
Loads a model's parameter dictionary (model.state_dict()) using a saved state_dict() object.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
torch.manual_seed(42)
epochs = 100 # Set the number of epochs
# Create empty loss lists to track values
train_loss_values = []
test_loss_values = []
epoch_count = []
for epoch in range(epochs):
### Training
model_0.train() # Put model in training mode (this is the default state of a model)
# 1. Forward pass on train data using the forward() method inside
y_pred = model_0(X_train)
# 2. Calculate the loss (how different are our models predictions to the ground truth)
loss = loss_fn(y_pred, y_train)
optimizer.zero_grad() # 3. Zero grad of the optimizer
loss.backward() # 4. Loss backwards
optimizer.step() # 5. Progress the optimizer
### Testing
# Put the model in evaluation mode
model_0.eval()
with torch.inference_mode():
# 1. Forward pass on test data
test_pred = model_0(X_test)
# 2. Caculate loss on test data
# predictions come in torch.float datatype, so comparisons need to be done with tensors of the same type
test_loss = loss_fn(test_pred, y_test.type(torch.float))
# Print out what's happening
if epoch % 10 == 0:
epoch_count.append(epoch)
train_loss_values.append(loss.detach().numpy())
test_loss_values.append(test_loss.detach().numpy())
print(f"Epoch: {epoch} | MAE Train Loss: {loss} | MAE Test Loss: {test_loss} ")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
torch.manual_seed(42)
# Set the number of epochs
epochs = 1000
# Put data on the available device
# Without this, error will happen (not all model/data on device)
X_train = X_train.to(device)
X_test = X_test.to(device)
y_train = y_train.to(device)
y_test = y_test.to(device)
for epoch in range(epochs):
### Training
model_1.train() # train mode is on by default after construction
# 1. Forward pass
y_pred = model_1(X_train)
# 2. Calculate loss
loss = loss_fn(y_pred, y_train)
# 3. Zero grad optimizer
optimizer.zero_grad()
# 4. Loss backward
loss.backward()
# 5. Step the optimizer
optimizer.step()
### Testing
model_1.eval() # put the model in evaluation mode for testing (inference)
# 1. Forward pass
with torch.inference_mode():
test_pred = model_1(X_test)
# 2. Calculate the loss
test_loss = loss_fn(test_pred, y_test)
if epoch % 100 == 0:
print(f"Epoch: {epoch} | Train loss: {loss} | Test loss: {test_loss}")
1 | class VAE(nn.Module): |
1 | # Create a neural net class |
1 | import torch |
1. https://github.com/mrdbourke/pytorch-deep-learning
2. https://readmedium.com/@rekalantar/variational-auto-encoder-vae-pytorch-tutorial-dce2d2fe0f5f
3. https://github.com/MicrosoftDocs/ml-basics/blob/master/05b%20-%20Convolutional%20Neural%20Networks%20(PyTorch).ipynb
4. https://github.com/ritchieng/the-incredible-pytorch?tab=readme-ov-file
5. https://github.com/claravania/lstm-pytorch
6. https://machinelearningmastery.com/pytorch-tutorial-develop-deep-learning-models/Suppose we have two sets of variable corresponding to two aspects such as height and weight, we want to analysis the relationship between this two sets. There are several ways to measure the relationship between them. However, sometime the it is hard to handle datasets with different dimensions, meaning, if \(X\in \mathbb{R}^m\) and \(Y\in \mathbb{R}^n\), how to resolve the relationship?
Assume there are two sets of data \(X\) and \(Y\), the size of \(X\) is \(n \times p\), whereas size of \(Y\) is \(n\times q\). That is, \(X\) and \(Y\) share the same row numbers but are differnt in columns number. The idea of CCA is simple: find the best match of \(X w_x\) and \(Y w_y\). Let's just set: \[X w_x = z_x\qquad\text{and}\qquad Y w_y = z_y\]
Where \(X\in \mathbb{R}^{n\times p}\), \(w_x \in \mathbb{R}^{p}\), \(z_x\in \mathbb{R}^n\), \(Y\in \mathbb{R}^{n\times q}\), \(w_y \in \mathbb{R}^{q}\), \(z_y\in \mathbb{R}^n\). \(w_x\) and \(w_y\) are often refered as canonical weight vectors, \(z_x\) and \(z_y\) are named images as well as canonical variates or canonical scores. To simplify the problem, we assume \(X\) and \(Y\) are standardized to zero mean and unit variance. Our task is to maximize the angle of \(z_x\) and \(z_y\), meaning:
\[\max_{z_x, z_y \in \mathbf{R^n}} <z_x, z_y>=\max \cos(z_x, z_y)=\max\frac{<z_x, z_y>}{\|z_x\|\|z_y\|}\]
with respect to: \(\|z_x\|_{2}=1\quad \|z_y\|_{2}=1\).
In fact, our task is just project \(X\) and \(Y\) to a new coordinate system after the linear transformation to \(X\) and \(Y\).
There are many solutions to this problems. Before start, We need make some assumptions: 1. the each column vector of \(X\) is perpendicular to the others. Which means \(X^T X= I\). The assumption is the same with \(Y\) and \(w_x, w_y\). We can find \(\min(p,q)\) canonical components, and the \(r\)th component is orthogonal to all the \(r-1\) components.
To solve the CCA problem using SVD, we first introduce the joint covariance matrix \(C\) such such that: \[\begin{equation} C = \begin{pmatrix} C_{xx} & C_{xy}\\ C_{yx} & C_{yy}\\ \end{pmatrix} \end{equation}\] Where \(C_{xx}=\frac{1}{n-1}X^\top X\) and \(C_{yy}=\frac{1}{n-1}Y^\top Y\) are the empirical variance matrices between \(X\) and \(Y\) respectively. The \(C_{xy}=\frac{1}{n-1} X^\top Y\) is the covariance matrix between \(X\) and \(Y\).
We next can reform CCA problem with two linear transformations \(w_x\) and \(w_y\) such that:
\[\begin{equation} w_x^\top C_{xx} w_x = I_p, \quad w_y^\top C_{yy} w_y = I_q, \quad w_x^\top C_{xy} w_y = D \end{equation}\] Where I_p and I_q are th p-dimensional and q-dimensional identity meatrics respectively. The diagonal matrix \(D = \text{diag}(\gamma_i)\) so that:
\[\begin{equation} \begin{pmatrix} {w}_x^\top & { 0}\\ { 0} & {w}_y^\top \end{pmatrix} \begin{pmatrix} C_{xx} & C_{xy}\\ C_{yx} & C_{yy} \end{pmatrix} \begin{pmatrix} {w}_x & { 0}\\ { 0} & {w}_y \end{pmatrix} = \begin{pmatrix} I_p & D\\ D^\top & I_q \end{pmatrix}, \end{equation}\]
The canoical variable: \[\begin{equation} Z_x = Xw_x, \quad Z_y = Y w_y \end{equation}\] The diagonal elements \(\gamma_i\) of D denote the canonical correlations. Thus we find the linear compounds \({Z}_x\) and \({Z}_y\) to maximize the cross-correlations. Since both \(C_{xx}\) and \(C_{yy}\) are symmetric positive definite, we can perform Cholesky Decomposition on them to get: \[\begin{equation} C_{xx} = C_{xx}^{\top/2} C_{xx}^{1/2}, \quad C_{yy} = C_{yy}^{\top/2} C_{yy}^{1/2} \end{equation}\]
where \(C_{xx}^{\top/2}\) is the transpose of \(C_{xx}^{1/2}\). Applying the inverses of the square root factors symmetrically on the joint covariance matrix \(C\), the matrix is transformed into: \[\begin{equation} \begin{pmatrix} C_{xx}^{-\top/2} & {\mathbf 0}\\ {\mathbf 0} & C_{yy}^{-\top/2} \end{pmatrix} \begin{pmatrix} C_{xx} & C_{ab}\\ C_{yx} & C_{yy} \end{pmatrix} \begin{pmatrix} C_{xx}^{-1/2} & {\mathbf 0}\\ {\mathbf 0} & C_{yy}^{-1/2} \end{pmatrix} = \begin{pmatrix} I_p & C_{xx}^{-1/2}C_{ab}C_{yy}^{-1/2}\\ C_{yy}^{-1/2}C_{yx}C_{xx}^{-1/2} & I_q \end{pmatrix}. \end{equation}\]
The canonical correlation problem is reduced to that of finding an SVD of a triple product: \[\begin{equation} U^{\top} (C_{xx}^{-1/2}C_{ab}C_{yy}^{-1/2}) V = D. \end{equation}\] The matrix \(C\) is thus reduced to the joint covariance matrix by applying a two-sided Jacobi method such that: \[\begin{equation} \begin{pmatrix} U^\top & {\mathbf 0}\\ {\mathbf 0} & V^\top \end{pmatrix} \begin{pmatrix} I_p & C_{xx}^{-1/2}C_{ab}C_{yy}^{-1/2}\\ C_{yy}^{-1/2}C_{_y}C_{xx}^{-1/2} & I_q \end{pmatrix} \begin{pmatrix} U & {\mathbf 0}\\ {\mathbf 0} & V \end{pmatrix} = \begin{pmatrix} I_p & D\\ D^\top & I_q \end{pmatrix} \end{equation}\]
with the desired transformation \({w}_x\) and \({w}_y\): \[\begin{equation} {w}_x = C_{xx}^{-1/2} U, \quad {w}_y = C_{yy}^{-1/2}V \end{equation}\] where the singular values \(\gamma_i\) are in descending order such that: \[\begin{equation} \gamma_1 \geq \gamma_2 \geq \cdots \geq 0. \end{equation}\]
The Problem can be reformed to solve the problem: \[\begin{equation} \underset{w_x \in \mathbb{R}^p, w_y\in \mathbb{R}^q}{\arg \max} w_x^\top C_{xy} w_y \end{equation}\] With respect to \(\|\|w_x^\top C_{xx} w_x\|\|_2 = \sqrt{w_x^\top C_{xx} w_x}=1\) and \(\|\|w_y^\top C_{yy} w_y\|\|_2 = \sqrt{w_y^\top C_{yy} w_y}=1\). The problem can apparently sovled by Lagrange multiplier technique. Let construct the Lagrange multiplier \(L\) such that: \[\begin{equation} L = w_x^\top C_{xy} w_y - \frac{\rho_1}{2} w_x^\top C_{xx} w_x - \frac{\rho_2}{2} w_y^\top C_{yy} w_y \end{equation}\]
The differentiation of L to \(w_x\) and \(w_y\) is: \[\begin{equation} \begin{aligned} \frac{\partial L}{\partial w_x} = C_{xy} w_y - \rho_1 C_{xx}w_x = \mathbf{0}\\ \frac{\partial L}{\partial w_y} = C_{yx} w_x - \rho_2 C_{yy}w_y = \mathbf{0} \end{aligned} \end{equation}\]
By left multipling \(w_x\) and \(w_y\) the above equation, we have:
\[\begin{equation} \begin{aligned} w_x^\top C_xy w_y -\rho_1 w_x^\top C_xx w_x = \mathbf{0}\\ w_y^\top C_yx w_x -\rho_2 w_y^\top C_yy w_y = \mathbf{0} \end{aligned} \end{equation}\] Since w_x^C_xx w_x = 1 and w_y^C_yy w_y = 1, we can obtain that \(\rho_1 = \rho_2 = \rho\). By substituting \(\rho\) to the formula. We can get: \[\begin{equation} w_x = \frac{C_{xx}^{-1}C_{xy}w_y}{rho} \end{equation}\] Evantually we have the equation: \[\begin{equation} C_{yx} C_{xx}^{-1} C_{xy} w_y = \rho^2 C_yy w_y \end{equation}\] Obviously, this is the form of eigenvalue decompostion problem where all eigen values are greater or equal to zero. By solving the eigenvalue decomposition we can find \(w_x\) and \(w_y\).
希望把博客的评论系统建立起来，之前使用的是disqus，重新部署的时候，页面大部分都无法显示。不想再用disqus。看到有人创造性的利用github作为载体建立评论系统，也就是Gitment了。
按照教程在github上设置了Gitment，惊闻Gitment需要请求服务，是作者搭的，作者已经不维护了。按照以下操作：
1
npm i --save gitment
修改自己js，连接自己搭建的服务器，WTF？
详细修改过程可参照：
https://sherry0429.github.io/2019/02/12/gitment%E4%BF%AE%E5%A4%8D/1
/node_modules/gitment/dist/gitment.browser.js
后继续寻觅其他可以评论系统，找到这篇文章： https://wangjiezhe.com/posts/2018-10-29-Hexo-NexT-3/ 根据此文章的教程安装了utterances。目前发现还是比较不错。知识现在看到的效果是全局评论。 issue-term不太了解具体，目前不想深入探究，仅仅设置pathname。
又一次，今年年初的又一次，博客系统hexo下的maupassant主题又罢工了。由于年初的各种事情繁琐而多，我就放弃治疗博客系统了，也就是说，有新的博文也无法发出来，先不管那些报错了。
现在稍微腾出一点儿时间，准备把博客系统好好弄一下。其实最简单的办法，也是屡试不爽的方法就是把所有的环境重新安装一遍，显示hexo，再是maupassant。这次不灵了，hexo
generate之后一堆报错。我甚至觉得maupassant已经无法搞定了，搜索错误的关键词，发现没有人遇到与我相同的问题。最后是怀疑我文章里有公式的特殊字符，影响markdown
parse。修改了hexo-renderer-marked的js，仍然有问题。最后决定换其他主题了，然而只要把所有的文章迁移过来一定会报错。报错如下：
1
2
3
4
5
6
7
8INFO Start processing
FATAL Something's wrong. Maybe you can find the solution here: http://hexo.io/docs/troubleshooting.html
Template render error: (unknown path) [Line 65, Column 565]
expected variable end
at Object._prettifyError (/Users/chengmingbo/blog_deploy/blog/node_modules/nunjucks/src/lib.js:36:11)
at Template.render (/Users/chengmingbo/blog_deploy/blog/node_modules/nunjucks/src/environment.js:542:21)
at Environment.renderString (/Users/chengmingbo/blog_deploy/blog/node_modules/nunjucks/src/environment.js:380:17
... ...
好吧，一切从头来，一个文件一个文件的添加，每次hexo generate一下。终于找到了一个有问题的文件。先注释掉再说，后面慢慢查是什么特殊字符引起的问题。好在可以更新博客了。
选定了 Next主题，又出现了反复，加上评论系统disqus发现博客白屏了，只有一个文件头显示。可是加功能一时爽，调试火葬场。当时实在记不起来到底是加了什么使得博客又不工作了。只能重头再来。在找主题的过程中发现star排名第四的hexo-theme-apollo已经停止开发，作者一句话让我决定不再折腾什么主题了：
1
专注文章内容的创作胜过博客样式的美观，祝各位玩的开心:
Pmf of Poisson Distribution is as follows:
\[f(X=k;\lambda)=\frac{\lambda^k e^{-\lambda}}{k!}\]
Our aim is to derive the the expectation of \(E(X)\) and the variance \(Var(X)\). Given that the formula of expectation: \[ E(X)=\sum_{k=0}^{\infty} k \frac{\lambda^k e^{-\lambda }}{k!} \]
Notice that when \(k=0\), the formula is equal to 0, that is:
\[\sum_{k=0}^{\infty} k \frac{\lambda^ke^{-\lambda}}{k!}\Large|_{k=0}=0\]
Then, the formula become as followed:
\[E(X)=\sum_{k=1}^{\infty} k \frac{\lambda^ke^{-\lambda}}{k!}\]
\[\begin{aligned}E(X)&=\sum_{k=0}^{\infty} k \frac{\lambda^ke^{-\lambda}}{k!}=\sum_{k=0}^{\infty} \frac{\lambda^ke^{-\lambda}}{(k-1)!}\\&=\sum_{k=0}^{\infty} \frac{\lambda^{k-1}\lambda e^{-\lambda}}{(k-1)!}\\&=\lambda e^{-\lambda}\sum_{k=1}^{\infty}\frac{\lambda^{k-1}}{(k-1)!}\end{aligned}\]
Now we need take advantage of Taylor Expansion, recall that:
\[e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots+\frac{x^{k-1}}{(k-1)!}=\sum_{k=1}^{\infty}\frac{x^{k-1}}{(k-1)!}\]
Compare \(E(X)\), we can get:
\[E(X)=\lambda e^{-\lambda}e^\lambda=\lambda\]
As known that \(Var(X)=E(X^2)-(E(x))^2\), we just get \(E(X^2)\). Given that:
\[E(X)=\sum_{k=1}^{\infty} k \frac{\lambda^ke^{-\lambda}}{k!}=\lambda\]
we can use this formula to derive the \(E(X^2)\),
\[\begin{aligned}E(X)=&\sum_{k=1}^{\infty} k \frac{\lambda^ke^{-\lambda}}{k!}=\lambda\\\Leftrightarrow&\sum_{k=1}^{\infty} k \frac{\lambda^k}{k!}=\lambda e^{\lambda}\\\Leftrightarrow&\frac{\partial\sum_{k=1}^{\infty} k \frac{\lambda^k}{k!}}{\partial \lambda}=\frac{\partial \lambda e^{\lambda}}{\partial \lambda}\\\Leftrightarrow&\sum_{k=1}^{\infty}k^2\frac{\lambda^{k-1}}{k!}=e^\lambda+\lambda e^\lambda\\\Leftrightarrow&\sum_{k=1}^{\infty}k^2\frac{\lambda^{k-1}e^{-\lambda}}{k!}=1+\lambda \\\Leftrightarrow&\sum_{k=1}^{\infty}k^2\frac{\lambda^{k}e^{-\lambda}}{k!}=\lambda+\lambda^2=E(X^2)\end{aligned}\]
then,
\[Var(X)=E(X^2)-(E(X))^2=\lambda+\lambda^2-(\lambda)^2=\lambda\]
Thus, we have proved that the Expectation and the Variance of Poisson Distribution are both \(\lambda\)
There are many classification algorithm such as Logistic Regression, SVM and Decision Tree etc. Today we'll talk about Gaussian Discriminant Analysis(GDA) Algorithm, which is not so popular. Actually, Logistic Regression performance better than GDA because it can fit any distributions from exponential family. However, we can learn more knowledge about gaussian distribution from the algorithm which is the most import distribution in statistics. Furthermore, if you want to understand Gaussian Mixture Model or Factor Analysis, GDA is a good start.
We, firstly, talk about Gaussian Distribution and Multivariate Gaussian Distribution, in which section, you'll see plots about Gaussian distributions with different parameters. Then we will learn GDA classification algorithm. We'll apply GDA to a dataset and see the consequnce of it.
As we known that the pdf(Probability Distribution Function) of
gaussian distribution is a bell-curve, which is decided by two
parameters \(\mu\) and \(\sigma^2\). The figure below shows us a
gaussian distribution with \(\mu=0\)
and \(\sigma^2=1\), which is often
referred to \(\mathcal{N}(\mu,
\sigma^2)\). Thus, Figure1 is distributed normally with \(\mathcal{N}(0,1)\). 
Figure 1. Gaussian Distribution with \(\mu=0\) and \(\sigma^2=1\).
Actually, parameter \(\mu\) and \(\sigma^2\) are exactly the mean and the variance of the distribution. Therefore, \(\sigma\) is the stand deviation of normal distribution. Let's take a look at area between red lines and magenta lines, which are respectively range from \(\mu\pm\sigma\) and from \(\mu\pm2\sigma\). The area between redlines accounts for 68.3% of the total area under the curve. That is, there are 68.3% samples are between \(\mu-\sigma\) and \(\mu+\sigma\) . Likely, there are 95.4% samples are between \(\mu-2\sigma\) and \(\mu+2\sigma\).
You must want to know how these two parameter influence the shape of PDF of gaussian distribution. First of all, when we change \(\mu\) with fixed \(\sigma^2\), the curve is the same as before but move along the random variable axis.
Figure 2. Probability Density Function curver with \(\mu=\pm2\) and \(\sigma=1\).
So, what if when we change \(\sigma\) then? Figure3. illustrates that smaller \(\sigma\) lead to sharper shape of pdf. Conversely, larger \(\sigma\) brings us broader curves.
Figure 3. Probability Density Function curver with \(\mu=0\) and change \(\sigma\).
Some may wonder what is the form of \(p(x)\) of a gaussian distribution, I just demonstrate here, you can compare Normal distribution with Multivariate Gaussian.
\[\mathcal{N(x|\mu, \sigma^2)}=\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2 } } \]
For convenience, we first see what is form of Multivariate Guassian Distribution:
\[\mathcal{N(x|\mu, \Sigma)}=\frac{1}{ { (2\pi)}^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } } e^{-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)}\]
where \(\mu\) is the mean, \(\Sigma\) is the covariance matrices, \(d\) is the dimension of random variable \(x\), specfically, 2-dimensional gaussian distribution, we have:
\[\mathcal{N(x|\mu, \Sigma)}=\frac{1}{\sqrt{2\pi}|\Sigma|^{\frac{1}{2 } } } e^{-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)}\]
In order to get an intuition of Multivariate Guassian Distribution, We first take a look at a distribution with \(\mu=\begin{pmatrix}0\\0\end{pmatrix}\) and \(\Sigma=\begin{pmatrix}1&0\\0&1\end{pmatrix}\).
Figure 4. 2-dimensional gaussian distribution with \(\mu=\begin{pmatrix}0\\0\end{pmatrix}\) and \(\Sigma=\begin{pmatrix}1&0\\0&1\end{pmatrix}\).
Notice that the the figure is rather than a curve but a 3-dimensional diagram. Just like normal distribution pdf, \(\sigma\) determines the shape of the figure. However, there are 4 entries of \(\Sigma\) can be changed in this example. Given that we need compute \(|\Sigma|\) as denominator and \(\Sigma^{-1}\) which demands non-zero determinant of \(\Sigma\), we must keep in mind that \(|\Sigma|\) is positive.
Rather than change \(\Sigma\), we firstly take a look at how the contour looks like when changing \(\mu\). Figure 5. illustrates the contour variation when changing \(\mu\). As we can see, we only move the center of the contour during the variation of \(\mu\). i.e. \(\mu\) detemines the position of pdf rather than the shape. Next, we will see how entries in \(\Sigma\) influence the shape of pdf.
Figure 5. Contours when change \(\mu\) with \(\Sigma=\begin{pmatrix}1&0\\0&1\end{pmatrix}\).
If scaling diagonal entries, we can see from figure 6. samples are concentrated to a smaller range when change \(\Sigma\) from \(\begin{pmatrix}1&0\\0&1\end{pmatrix}\) to \(\begin{pmatrix}0.3&0\\0&0.3\end{pmatrix}\). Similarly, if we alter \(\Sigma\) to \(\begin{pmatrix}3&0\\0&3\end{pmatrix}\), then figure will spread out.
Figure 6. Density when scaling diagonal entries to 0.3.
What if we change only one entry of the diagonal? Figure 7. shows the variation of the density when change \(\Sigma\) to \(\begin{pmatrix}1&0\\0&5\end{pmatrix}\). Notice the parameter spuashes and stretches the figure along coordinate axis.
Figure 7. Density when scaling one of the diagonal entries.
We now try to change entries along secondary diagonal. Figure 8. demonstrates that the variation of density is no longer parallel to \(X\) and \(Y\) axis, where \(\Sigma=\begin{pmatrix}1 &0.5\\0.5&1\end{pmatrix}\).
Figure 8. Density when scaling secondary diagonal entries to 0.5
When we alter secondary entries to negative 0.5, the direction of contour presents a mirror to contour when positive.
Figure 9. Density when scaling secondary diagonal entries to -0.5
In light of the importance of determinant of \(\Sigma\), what will happen if the determinant is close to zero. Actually, we can, informally, take determinant of a matrice as the volume of which. Similarly, when determinant is smaller, the volume under density curve become smaller. Figure 10. illustrates the circumstance we talked above where \(\Sigma=\begin{pmatrix}1&0.99\\0.99&1\end{pmatrix}\).
Figure 10. Density when determinant is close to zero.
When input features \(x\) are continuous variables, we can use GDA classify data. Firstly, let's take a look at how GDA to do the job. Figure 11. show us two gaussian distributions, they share the same covariance \(\Sigma=\begin{pmatrix}1&0\\0&1\end{pmatrix}\) , and repectively with parameter \(\mu_0=\begin{pmatrix}1\\1\end{pmatrix}\) and \(\mu_1=\begin{pmatrix}-1\\-1\end{pmatrix}\). Imagine you have some data which fall into the cover of the first and second Gaussian Distribution. If we can find such distributions to fit the data, then we'll have the capcity to decide which is new data coming from, the first or the second one.
Figure 11. Two gaussian distributions with respect to \(\mu_0=\begin{pmatrix}1\\1\end{pmatrix}\) and \(\mu_1=\begin{pmatrix}-1\\-1\end{pmatrix}\) , and \(\Sigma=\begin{pmatrix}1&0\\0&1\end{pmatrix}\)
Specifically, let's look at a concrete example, Figure 12 are samples drawn from two Gaussian distribution. There are 100 blue '+'s and 100 red 'o's. Assume that we have such data to be classified. We can apply GDA to solve the problem.
CODE:
1 | %% octave |
Figure 12. 200 samples drawn from two Gaussian Distribution with parameters \(\mu_0=\begin{bmatrix}1\\1\end{bmatrix},\mu_1=\begin{bmatrix}4\\4\end{bmatrix},\Sigma=\begin{bmatrix}1&0\\0&1\end{bmatrix}\).
Now, let's define the algorithm. Firstly we assume discrete random variable classes \(y\) are distributed Bernoulli and parameterized by \(\phi\), then we have:
\[y\sim {\rm Bernoulli}(\phi)\]
Concretely, the probablity of \(y=1\) is \(\phi\), and \(1-\phi\) when \(y=0\). We can simplify two equations to one:
\(p(y|\phi)=\phi^y(1-\phi)^{1-y}\)
Apparently, \(p(y=1|\phi)=\phi\) and \(p(y=0|\phi)=1-\phi\) given that y can only be \(0\) or \(1\).
Another assumption is that we consider \(x\) are subject to different Gaussian Distributions given different \(y\). We assume the two Gaussian distributions share the same covariance and different \(\mu\). Based on above all, then
\[p(x|y=0)=\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x-\mu_0)^{T}\Sigma^{-1}(x-\mu_0)}\]
\[p(x|y=1)=\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x-\mu_1)^{T}\Sigma^{-1}(x-\mu_1)}\]
i.e. \(x|y=0 \sim \mathcal{N}(\mu_0,\Sigma)\) and \(x|y=1 \sim \mathcal{N}(\mu_1,\Sigma)\). suppose we have \(m\) samples, it is hard to compute \(p(x^{(1)}, x^{(2)}, x^{(3)},\cdots,x^{(m)}|y=0)\) or \(p(x^{(1)}, x^{(2)}, x^{(3)},\cdots,x^{(m)}|y=1)\) . In general, we assume the probabilty of \(x^{(i)}\) \(p(x^{(i)}|y=0)\) is independent to any \(p(x^{(j)}|y=0)\), then we have:
\[p(X|y=0)=\prod_{i=1\,y^{(i)}\neq1}^{m}\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x^{(i)}-\mu_0)^{T}\Sigma^{-1}(x^{(i)}-\mu_0)}\]
Vice versa,
\[p(X|y=1)=\prod_{i=1\,y^{(i)}\neq 0}^{m}\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x^{(i)}-\mu_1)^{T}\Sigma^{-1}(x^{(i)}-\mu_1)}\]
Here \(X=(x^{(1)}, x^{(2)}, x^{(3)},\cdots,x^{(m)})\). Now, we want to maximize \(p(X|y=0)\) and \(p(X|y=1)\). Why is that, because we hope find parameters that let \(p(X|y=0)p(X|y=1)\) largest, based on that the samples are from the two Gaussian Distributions. These samples we have are more likely emerging. Thus, our task is to maximize \(p(X|y=0)p(X|y=1)\) , we let
\[\mathcal{L}(\phi,\mu_0,\mu_1,\Sigma)=\arg\max p(X|y=0)p(X|y=1)=\arg\max\prod_{i=1}^{m}p(x^{(i)}, y^{(i)};\phi,\mu_0,\mu_1,\Sigma)\]
It's tough for us to maximize \(\mathcal{L}(\phi,\mu_0,\mu_1,\Sigma)\). Notice function \(\log\) is monotonic increasing. Thus, we can maximize \(\log\mathcal{L}(\phi,\mu_0,\mu_1,\Sigma)\) instead of \(\mathcal{L}(\phi,\mu_0,\mu_1,\phi)\), then:
\[\begin{aligned}\ell(\phi,\mu_0,\mu_1,\Sigma)&=\log\mathcal{L}(\phi,\mu_0,\mu_1,\Sigma)\\&=\arg\max\log\prod_{i=1}^{m}p(x^{(i)},y^{(i)};\phi,\mu_0,\mu_1,\Sigma)\\&=\arg\max\log\prod_{i=1}^{m}p(x^{(i)}|y^{(i)};\mu_0,\mu_1,\Sigma)p(y^{(i)};\phi)\\&=\arg\max\sum_{i=1}^{m}p(x^{(i)}|y^{(i)};\mu_0,\mu_1,\Sigma)+p(y^{(i)};\phi)\\&=\arg\max\sum_{i=1}^{m}p(x^{(i)}|y^{(i)};\mu_0,\mu_1,\Sigma)+\sum_{i=1}^{m}p(y^{(i)};\phi)\end{aligned}\]
By now, we have found a convex function with respect parameters \(\mu_0, mu_1,\Sigma\) and \(\phi\). Next section, we'll obtain these parameter through partial derivative.
To estimate these four parameters, we just apply partial derivative to \(\ell\). Now we estimate \(\phi\) in the first place. We let \(\frac{\partial \ell}{\partial \phi}=0\), then
\[\begin{aligned}\frac{\partial \ell(\phi,\mu_0,\mu_1,\Sigma)}{\partial \phi}=0&\Rightarrow\frac{\partial \arg\max\sum_{i=1}^{m}p(x_i|y;\mu_0,\mu_1,\Sigma)+\sum_{i=1}^{m}p(y_i;\phi)}{\partial \phi}=0\\&\Rightarrow\frac{\partial\sum_{i=1}^{m}\log p(y^{(i)};\phi)}{\partial \phi}=0\\&\Rightarrow\frac{\partial\sum_{i=1}^{m}\log \phi^{y^{(i) } } (1-\phi)^{(1-y^{(i)}) } } {\partial \phi}=0\\&\Rightarrow\frac{\partial\sum_{i=1}^{m}{y^{(i) } } \log \phi+{(1-y^{(i)})}\log(1-\phi)}{\partial \phi}=0\\&\Rightarrow\frac{\partial\sum_{i=1}^{m}{ { 1}{\{y^{(i)}=1\} } }\log \phi+{1}{\{y^{(i)}=0\} } \log(1-\phi)}{\partial \phi}=0\\&\Rightarrow\phi=\frac{1}{m}\sum_{i=1}^{m}1\{y^{(i)}=1\}\end{aligned}\]
Note that \(\mu_0\) and \(\mu_1\) is symmetry in the equation, thus, we need only obtain one of them. Here we take the derivative to \(\mu_0\)
\[\begin{aligned}\frac{\partial \ell(\phi,\mu_0,\mu_1,\Sigma)}{\partial \mu_0}=0&\Rightarrow\frac{\partial \arg\max\sum_{i=1}^{m}p(x_i|y;\mu_0,\mu_1,\Sigma)+\sum_{i=1}^{m}p(y_i;\phi)}{\partial \mu_0}=0\\&\Rightarrow\frac{\partial\sum_{i=1}^{m} \log p(x^{(i)}|y^{(i)};\mu_0,\mu_1,\Sigma)}{\partial \mu_0}=0\\&\Rightarrow\frac{\partial \sum_{i=1}^{m}\log\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x^{(i)}-\mu_0)^T\Sigma^{-1}(x^{(i)}-\mu_0) } } {\partial \mu_0}=0\\&\Rightarrow0+\frac{\partial \sum_{i=1}^{m}{-\frac{1}{2}(x^{(i)}-\mu_0)^T\Sigma^{-1}(x^{(i)}-\mu_0) } } {\partial \mu_0}=0\end{aligned}\]
We have \(\frac{\partial X^TAX}{\partial X}=(A+A^T)X\)，let \((x^{(i)}-\mu_0)=X\), then,
\[\begin{aligned}\frac{\partial \ell(\phi,\mu_0,\mu_1,\Sigma)}{\partial \mu_0}=0&\Rightarrow 0+\frac{\partial \sum_{i=1}^{m}{-\frac{1}{2}(x^{(i)}-\mu_0)^T\Sigma^{-1}(x^{(i)}-\mu_0) } } {\partial \mu_0}=0\\&\Rightarrow{\sum_{i=1}^{m}-\frac{1}{2}((\Sigma^{-1})^T+\Sigma^{-1})(x^{(i)}-\mu_0)\cdot(-1)}=0\\&\Rightarrow \sum_{i=1}^{m}1\{y^{(i)}=0\}x^{(i)}=\sum_{i=1}^{m}1\{y^{(i)}=0\}\mu_0\\&\Rightarrow\mu_0=\frac{\sum_{i=1}^{m}1\{y^{(i)}=0\}x^{(i) } } {\sum_{i=1}^{m}1\{y^{(i)}=0\} } \end{aligned}\]
Simlarly,
\[\mu_1=\frac{\sum_{i=1}^{m}1\{y^{(i)}=1\}x^{(i) } } {\sum_{i=1}^{m}1\{y^{(i)}=1\} } \]
Before calculate \(\Sigma\), I first illustrate the truth that \(\frac{\partial|\Sigma|}{\partial\Sigma}=|\Sigma|\Sigma^{-1},\quad \frac{\partial\Sigma^{-1 } } {\partial\Sigma}=-\Sigma^{-2}\), then
\[\begin{aligned}\frac{\partial \ell(\phi,\mu_0,\mu_1,\Sigma)}{\partial \Sigma}=0&\Rightarrow\frac{\partial\sum_{i=1}^{m} \log p(x^{(i)}|y^{(i)};\mu_0,\mu_1,\Sigma)+\sum_{i=1}^{m} \log p(y^{(i)};\phi)}{\partial \Sigma}=0\\&\Rightarrow\frac{\partial \sum_{i=1}^{m}\log\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}(x^{(i)}-\mu_{y^{(i) } } )^T\Sigma^{-1}(x^{(i)}-\mu_{y^{(i) } } ) } } {\partial \Sigma}=0\\&\Rightarrow\frac{\partial \sum_{i=1}^{m}-\frac{d}{2}\log2\pi}{\partial \Sigma}+\frac{\partial \sum_{i=1}^{m}-\frac{1}{2}\log|\Sigma|}{\partial \Sigma}+\frac{\partial \sum_{i=1}^{m}{-\frac{1}{2}(x^{(i)}-\mu_{y^{(i) } } )^T\Sigma^{-1}(x^{(i)}-\mu_{y^{(i) } } ) } } {\partial \Sigma}=0\\&\Rightarrow\frac{\partial \sum_{i=1}^{m}-\frac{1}{2}\log|\Sigma|}{\partial \Sigma}+\frac{\partial \sum_{i=1}^{m}{-\frac{1}{2}(x^{(i)}-\mu_{y^{(i) } } )^T\Sigma^{-1}(x^{(i)}-\mu_{y^{(i) } } ) } } {\partial \Sigma}=0\\&\Rightarrow m\frac{1}{|\Sigma|}|\Sigma|\Sigma^{-1}+\sum_{i=1}^m(x^{(i)}-\mu_{y^{(i) } } )^T(x^{(i)}-\mu_{y^{(i) } } )(-\Sigma^{-2}))=0\\&\Rightarrow\Sigma=\frac{1}{m}\sum_{i=1}^{m}(x^{(i)}-\mu_{y^{(i) } } )(x^{(i)}-\mu_{y^{(i) } } )^T\end{aligned}\]
In spite of the harshness of the deducing, the outcome are pretty beautiful. Next, we will apply these parameters and see how the estimation performance.
Notice the data drawn from two Gaussian Distribution is random, thus, if you run the code, the outcome may be different. However, in most cases, distributions drawn by estimated parameters are roughly the same as the original distributions.
\[\begin{aligned}&\phi=0.5\\&\mu_0=\begin{bmatrix}4.0551\\4.1008\end{bmatrix}\\&\mu_1=\begin{bmatrix}0.85439\\1.03622\end{bmatrix}\\&\Sigma=\begin{bmatrix}1.118822&-0.058976\\-0.058976&1.023049\end{bmatrix}\end{aligned}\]
From Figure 13, We can see contours of two Gaussian distribution, and most of samples are correctly classified.
CODE:
1 | %% octave |
Figure 13. Contours drawn from parameters estimated.
In fact, we can compute the probability of each data point to predict which distribution it is more likely belongs, for example, if we want to predict \(x=\begin{pmatrix}0.88007\\3.9501\end{pmatrix}\) is more of the left distribution or the right, we apply \(x\) to these two distribution:
\[\begin{aligned}p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg|y=0\right)=&\frac{1}{2\pi|\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}{\begin{bmatrix}x_1-\mu_1\\x_2-\mu_2\end{bmatrix}^T\Sigma^{-1}\begin{bmatrix}x_1-\mu_1\\x_2-\mu_2\end{bmatrix} } }\\=&\frac{1}{ { 2\pi}\left|\begin{matrix}1.1188&-0.059\\-0.059&1.023\end{matrix}\right|^{\frac{1}{2 } } }e^{-\frac{1}{2}{\begin{bmatrix}-3.175\\-0.151\end{bmatrix}^T\begin{bmatrix}0.896&-0.052\\-0.0520&0.98\end{bmatrix}\begin{bmatrix}-3.175\\-0.151\end{bmatrix} } }\\=&\frac{1}{2\pi\sqrt{(1.141) } } e^{-\frac{1}{2}\times 9.11}=0.149\times 0.01=0.0015\end{aligned}\]
and
\[\begin{aligned}p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg| y=1\right)&\frac{1}{2\pi|\Sigma|^{\frac{1}{2 } } }e^{-\frac{1}{2}{\begin{bmatrix}x_1-\mu_1\\x_2-\mu_2\end{bmatrix}^T\Sigma^{-1}\begin{bmatrix}x_1-\mu_1\\x_2-\mu_2\end{bmatrix} } }\\=&\frac{1}{ { 2\pi}\left|\begin{matrix}1.1188&-0.059\\-0.059&1.023\end{matrix}\right|^{\frac{1}{2 } } }e^{-\frac{1}{2}{\begin{bmatrix}0.03\\2.91\end{bmatrix}^T\begin{bmatrix}0.896&-0.052\\-0.0520&0.98\end{bmatrix}\begin{bmatrix}0.03\\2.91\end{bmatrix} } }\\=&\frac{1}{2\pi\sqrt{(1.141) } } e^{-\frac{1}{2}\times 8.336}=0.149\times 0.015=0.0022\end{aligned}\]
In light of the equivalency of \(p(y=1)\) and \(p(y=0)\) (both are \(0.5\)), we just compare\(p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg| y=1\right)\) to \(p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg| y=0\right)\). Apparently, this data point is predicted from the left distribution, which is a wrong assertion. Actually, in this example, we have only this data pointed classified incorrectly.
You may wonder why there is a blue line. It turns out that all the data point below the blue line will be considered as blue class. Otherwise, data points above the line is classified as the red class. How it work?
The blue line is decision boundary, if we know the expression of this line, the decision will be made easier. In fact GDA is a linear classifier, we will prove it later. Still, we see the data point above, if we just divide one probability to another, we just need find if the ratio larger or less than 1. For our example, the ratio is roughly 0.68, so the data point is classified to be the blue class.
\[\frac{p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg| y=0\right)p(y=0)}{p\left(x=\begin{bmatrix}0.88\\3.95\end{bmatrix}\Bigg|y=1\right)p(y=1)}=\frac{0.0015}{0.0022}=0.68182<1\]
Figure 14. Decision Boundary
If we can obtain the expression of the ratio, that should be good. So given a new \(x\), we predict problem is tranformed as followed:
\[x\in \text{red class}\propto\mathcal{R}=\frac{p(x|y=1)p(y=1)}{p(x|y=0)p(y=0)} > 1\]
which is equal to
\[\mathcal{R}=\log\frac{p(x|y=1)p(y=1)}{p(x|y=0)p(y=0)} =\log\frac{\phi}{1-\phi}+\log\frac{\mathcal{N}(x;\mu_1,\Sigma)}{\mathcal{N}(x;\mu_0,\Sigma)}> 0\]
Then,
\[\begin{aligned}\mathcal{R}&=\log\frac{\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }\exp(-\frac{1}{2}(x-\mu_1)^T\Sigma^{-1}(x-\mu_1))}{\frac{1}{(2\pi)^{\frac{d}{2 } } |\Sigma|^{\frac{1}{2 } } }\exp(-\frac{1}{2}(x-\mu_0)^T\Sigma^{-1}(x-\mu_0))}+\log\frac{\phi}{1-\phi}\\&=-\frac{1}{2}(x-\mu_1)^T\Sigma^{-1}(x-\mu_1))+\frac{1}{2}(x-\mu_0)^T\Sigma^{-1}(x-\mu_0))+\log\frac{\phi}{1-\phi}\\&=-\frac{1}{2}x^T\Sigma^{-1}x+\mu_1^T\Sigma^{-1}x-\frac{1}{2}\mu_1^T\Sigma^{-1}\mu_1+-\frac{1}{2}x^T\Sigma^{-1}x-\mu_0^T\Sigma^{-1}x+\frac{1}{2}\mu_0^T\Sigma^{-1}\mu_0+\log\frac{\phi}{1-\phi}\\&=(\mu_0-\mu_1)^T\Sigma^{-1}x-\frac{1}{2}\mu_1^T\Sigma^{-1}\mu_1+\frac{1}{2}\mu_0^T\Sigma^{-1}\mu_0+\log\frac{\phi}{1-\phi}\end{aligned}\]
Here, \(\mu_1^T\Sigma^{-1}x=x^T\Sigma^{-1}\mu_1\) because it is a real number. For a real number \(a=a^T\), moreover, \(\Sigma^{-1}\) is symmetric, so \(\Sigma^{-T}=\Sigma^{-1}\). Let's set \(w^T=(\mu_1-\mu_0)^T\Sigma^{-1}\) and \(w_0=-\frac{1}{2}\mu_1^T\Sigma^{-1}\mu_1+\frac{1}{2}\mu_0^T\Sigma^{-1}\mu_0+\log\frac{\phi}{1-\phi}\), then we have:
\[\mathcal{R}=\log\frac{p(x|y=1)p(y=1)}{p(x|y=0)p(y=0)} =w^Tx+w_0\]
If you plug parameters in the formula, you will find:
\(\mathcal{R}=-3.0279x_1-3.1701x_2+15.575=0\)
It is the decision boundary(Figure 14.). Since you have got the decision boundary formula, it is convenient to use the decision boundary function predict if a data point \(x\) belongs to the blue or red class. If \(\mathcal{R}>0\), \(x\in \text{red class}\), otherwise, \(x\in \text{blue class}\).
Today, we have talked about Guassian Distribution and its Multivariate form. Then, we assume two groups of data drawn from Gaussian Distributions. We apply Gaussian Discriminant Analysis to the data. There are 200 data point, only one is misclassified. In fact we can deduce GDA to Logistic regression Algorithm(LR). But LR can not deduce GDA, i.e. LR is a better classifier, especially when we do not know the distribution of the data. However, if you have known that data is drawn from Gaussian Distribution, GDA is the better choice.