*Andrea Arauza Rivera, Ph. D.

Assistant Professor, Mathematics

Cal State East Bay

Matt McClinton

Cal State East Bay

David Smith

Cal State East Bay

*Corresponding author.

Numberlink is a puzzle game in which players are given a grid with nodes marked with a natural number, $n$, and asked to create $n$ connections with neighboring nodes. Connections can only be made with top, bottom, left and right neighbors, and one cannot have more than two connections between any neighboring nodes.
In this paper, we give a mathematical formulation of the puzzles via graphs and give some immediate consequences of this formulation.

The main result of this work is an algorithm which provides insight into characteristics of these puzzles and their solutions.

Finally, we give a few open questions and further directions.

§ INTRODUCTION

With millions of downloads between the Apple and Google Play stores, Puzzledom is an app that offers a variety of curious puzzles. One of these is the puzzle called Numberlink [1]. Puzzles like Numberlink have been the subject of a number of interesting articles; the reader may enjoy any one of the following [2], [3]. The reader should note that there is another popular game called Numberlink.
This version of the game is discussed in [4], [5], [6].

We focus on the Numberlink puzzles found in the Puzzledom app. These puzzles begin with a set of numbered boxes (nodes) configured in a grid; see Figure <ref>. The puzzle is solved when the player creates links between the numbered nodes so that

* a node with number $n$ has $n$ connections to other nodes,

\[ \includegraphics[scale = .35]{pics/miniex1}\]

* there are no more than $k=2$ connections between any two nodes, and

\[ \includegraphics[scale = .35]{pics/miniex2}\]

* the connections must create a path between any two given nodes (path connected).

\[ \includegraphics[scale = .35]{pics/miniex3}\]

* The game also implicitly requires that connections be made horizontally or vertically, and that no connections intersect.

\[ \includegraphics[scale = .35]{pics/miniex4}\]

Three sample Numberlink puzzles are shown in Figure <ref>. The reader is encouraged to whip out a pencil and try to solve each puzzle!

We show 3 sample puzzles for the reader to try and enjoy. Each of these puzzles is included in Numberlink and listed as “Novice” (left) or “Regular”(center and right) [1].

The sections in this article are set up as follows:

* Section 1 gives a mathematical formulation of the Numberlink game via graphs. This section also includes some initial consequences of this formulation.

* Section 2 describes ways in which the player may find the best nodes to start to create connections. This analysis is based on the number in the node as well as the number of neighbors available for connecting.

* Section 3 contains the main results of this article; an algorithm which produces the “guaranteed connections” between nodes. In this section we prove that if the algorithm we outline arrives at a solution, then the solution is unique.

* We conclude in section 4 with some closing thoughts and open questions about these puzzles.

§ NUMBERED $K$-GRIDS

This section contains the definitions needed to describe the Numberlink puzzles in terms of graphs.  We begin by defining the initial set up of a puzzle as a graph where the nodes are labeled with a whole number $n$ and arranged on a grid. We call these numbered $k$-grids. Next, we define what it means for a node $p$ in a numbered $k$-grid to have top, bottom, left and right neighbors. Finally, we define what it means for a numbered $k$-grid to be solved.

The reader may be wondering what on earth this $k$ business is. Indeed, in the original Puzzledom-Numberlink puzzles there is a rule that no two nodes may share more than $k=2$ connections. We work with a more general rule and allow $k$ to be any positive whole number.

A numbered k-grid is a finite collection of nodes and connections between nodes satisfying the following:

* each node is given coordinates $(x, y)$ where $x, y \in \N \cup \{0\}$;

* each node is labeled with a
magnitude $n\in \N$;

* there are no more than $k$ connections between nodes;

* connections can only exist connecting horizontally or vertically adjacent nodes.

We denote a numbered $k$-grid by $\Gamma_k$ and write $p = [x, y, n]$ for a node in $\Gamma_k$. Often we will denote the magnitude of $p$ as $magn(p)$.

A typical initial set-up for a Numberlink game is shown (left). Two examples of fully connected nodes, given this initial set-up, are shown (center-left and center-right). Finally, a full solution of the game with the initial set-up is shown (right).

The choice of a square to represent the nodes in a numbered grid is somewhat arbitrary. One could use a circle or other figure to represent nodes. We choose to use a square to emphasize the top, bottom, left and right neighbors of the node.
We now define what it means to be neighboring to a node $p$.

Neighboring Nodes in $\Gamma_k$

Let $W = \{1,2,3,4\}$ and $\Gamma_k$ be a numbered $k$-grid. Consider any node $p = [x_p,y_p,n_p] \in \Gamma_k$. We refer to $1_p,2_p,3_p,4_p$ as the top, right, bottom, and left neighbors to p, respectively. More formally, for each $l \in W$, we define the $l_p$ neighbor by the following:

\begin{align*}
&1_p = [x_{1_p}, \ y_{1_p}, \ n_{1_p}] \in \Gamma_k \text{ where } x_{1_p}=x_{p}, \text{ and } y_{1_p} = \min \{y : [x_{p}, y, m] \in \Gamma_k \text{ and }y > y_p \}. \\
& \newline \\
&2_p = [x_{2_p}, \ y_{2_p}, \ n_{2_p}] \in \Gamma_k \text{ where } y_{2_p}=y_{p}, \text{ and }x_{2_p} = \min \{x : [x, y_p, m] \in \Gamma_k \text{ and }x > x_p \}. \\
& \newline \\
&3_p = [x_{3_p}, \ y_{3_p}, \ n_{3_p}] \in \Gamma_k \text{ where } x_{3_p}=x_{p}, \text{ and } y_{3_p} = \max \{y : [x_p, y, m] \in \Gamma_k \text{ and }y < y_p \}. \\
& \newline \\
&4_p = [x_{4_p}, \ y_{4_p}, \ n_{4_p}] \in \Gamma_k \text{ where } y_{4_p}=y_{p}, \text{ and } x_{4_p} = \max \{x : [x, y_p, m] \in \Gamma_k \text{ and }x < x_p \}.
\end{align*}

If no node in $\Gamma_k$ satisfies the conditions of an $l_p$ neighbor, we say the $l_p$ neighbor does not exist.

In Numberlink, when one creates a connection between nodes, the magnitude of each of the two nodes is reduced by 1. A puzzle is completed when the magnitude of all nodes has been reduced to 0 and no rules of the puzzle have been violated (e.g. no nodes share more than $k$ connections, the graph is path connected, connections are only shared by neighboring nodes, and no connections cross).
We now give the definition of a solved $k$-grid in our context of graphs.

A solved $k$-grid

For any numbered $k$-grid, $\Gamma_k$, we say it is solved if the set of connections in $\Gamma_k$ satisfy the following:

* Every node $p \in \Gamma_k$ is a completed node, meaning the number of connections between $p$ and it's neighboring nodes equals $magn(p)$.

* The number of connections between any two neighboring nodes is at most $k$.

* No pair of connections

* The graph $\Gamma_k$ is path connected.

§.§ Some initial results

Three examples of the ways in which $k$ affects the number of solutions that can exist for a $k$-grid.

Every puzzle in Numberlink can be classified as a solvable $2$-grid.
Note also that any solvable $1$-grid is also a solvable $2$-grid. More generally, if $\Gamma_k$ is solvable, then the graph with the same nodes as $\Gamma_k$ can be solved with any $k'>k$.

Let's take a moment to inspect how $k$ influences solvability. Consider the 3 examples at the top of Figure <ref>.

We note the following:

* For $k=1$, only the left grid is a solvable $1$-grid.

* Moving to $k=2$ with the left grid,
the graph resulting from making 2 connections between nodes would be disjoint. Hence the left grid is a solvable $2$-grid if and only if a single connection is formed between all nodes.

* With $k=2$, the middle and right grids are now solvable. The middle grid has two solutions and the right grid is uniquely solved.

* If $k=3$, we see that the right grid has another 2 solutions. No new solutions are available for the left or middle grids.

The above are examples of $k$-grids whose solutions are relatively simple to build by inspection. As one progresses through the vast array of puzzles in Numberlink, eyeballing solutions involves a good deal of mental planning. Every puzzle in Numberlink has at least one solution, yet an arbitrary $k$-grid has no guarantee of being solvable. Are there ways of quickly identifying unsolvable $k$-grids by inspection? In fact, there are!

Let $\Gamma_k$ be any numbered $k$-grid. If any of the following hold, then $\Gamma_k$ is unsolvable.

* $\Gamma_k$ has a node with no neighbors.

* The sum of the magnitudes of all nodes in $\Gamma_k$ is odd.

* There is a node $p \in \Gamma_k$ such that the sum of the magnitudes of the neighbors of $p$ is less than the magnitude of $p$.

* Every configuration of connections for some node $p$ results in a disjoint graph.

* A node $p$ has $r$ neighbors $(\text{where } r = 1,2,3,4)$ and $magn(p) > rk$.

* (Incompatible nodes property) If $k>1$ and there exists a node $p$ with $r$ neighboring nodes such that $magn(p)=(r-1)k+j$ for $j=2,3,\dots,k$, and $p$ has a neighboring node $q$ with $magn(q) \leq j-1$.

For part 1, if a node in any $k$-grid has no neighboring nodes, then it is impossible to form any connections, and thus unsolvable for any $k$.

Part 2 holds true by the “first theorem of graph theory” [7]
. Since each connection reduces the magnitudes of the two corresponding nodes by 1, there will be no way to complete every node while also having the sum of the magnitudes be odd.

For part 3, let $\Gamma_k$ be some $k$-grid with a node $p$ where the sum of the magnitudes of the neighboring nodes is less than the magnitude of $p$, i.e
$$magn(p) > \sum\limits_{l \in W} magn(l_p).$$

To say $\Gamma_k$ is solvable means there exist a connection configuration for $p$ where $p$ is complete. Yet the neighboring nodes to $p$ cannot take in $magn(p)$ connections.

The proof for part 4 also follows from the definition of a solved $k$-grid. Let $\Gamma_k$ be a $k$-grid with a node $p$ such that every edge formation results in $\Gamma_k$ disjoint. Meaning all completions of $p$ violates the criteria of being solved and hence $\Gamma_k$ unsolvable.

Now consider any $k$-grid $\Gamma_k$ satisfying the condition in part 5 for $r=1$. To say $p$ has one neighbor with $magn(p) > k$ implies $magn(p) \geq k+1$. Meaning that there must be more than $k$ connections between $p$ and its neighboring node in $\Gamma_k$, and hence $\Gamma_k$ is unsolvable. For $r=2$, $magn(p) > 2k$ implies $magn(p) \geq 2k+1$. Again this means that more than $k$ connections must be formed for at least one neighboring node, which deems $\Gamma_k$ as unsolvable. Similar arguments can be made for $r=3,4$.

For part 7, if $magn(p) = (r-1)k+j$ with $j = 2,3,4, \dots,k$, then $p$ will need to have $k$ connections with $(r-1)$ of it's neighbors and $j$ connections with the remaining neighbor. Thus if $p$ has a neighbor $q$ with $magn(q) \leq j-1$, then there will be no way to complete the node $p$.
$\quad \blacksquare$

The reader is encouraged to find their own condition which guarantees that a $k$-grid will be unsolvable. Note that the incompatible nodes property tells us about what magnitudes two neighboring nodes cannot have and still allow the grid to be solvable. For example, in a 2-grid, a node of magnitude 8 with four neighbors cannot have a neighbor of magnitude 1.
The reader should take a moment to write out other instances in which the incompatible nodes property tells us that a $k$-grid will be unsolvable.

§.§ Creating connections between nodes

Now let's talk about connections. Nodes of relatively high and low magnitudes often have little variety in how the connections are constructed. A node of magnitude 3 with four neighbors in a $2$-grid has only 4 ways of creating connections. However, a node of magnitude 20 with four neighbors in a $10$-grid has 891 ways of drawing connections! It is not always the case that all neighboring nodes can take in a full $k$ connections. This means that for the magnitude 20 node, we don't always have to consider all 891 cases. So how do we avoid exhausting every feasible case? We will now begin to describe an algorithm to determine what connections a node must have.

Connections between Neighboring Nodes

Let $\Gamma_k$ be a $k$-grid and $p \in \Gamma_k$. Let $W = \{1,2,3,4\}$. Define $\Phi_k[p]$ to be the set of all connection configurations for the node $p$ within $\Gamma_k$.

Thus we have

$$\Phi_k[p] = \{ \omega = \omega_1 \omega_2 ... \omega_n : \omega_i \in W, n=magn(p), \omega_i = \omega_j \text{ at most k times}\}. $$

Here the symbols $1, 2, 3, 4$  refer to the top, right, bottom or left neighbor of the node $p$, respectively.
We denote the length of a word by $||\omega||$. We say two words are equivalent, $\omega = \omega'$ in $\Phi_k[p]$, when they have the same symbols without deleting duplicates $\{\omega_i\} = \{\omega'_j\}$.

We emphasize that the rearrangement of the $\omega_i$ does not result in a new word. To say $\omega = 1212121$ means that 4 connections are created with the top ($\omega_i = 1$) neighbor, and 3 connections are created with the right ($\omega_i = 2$) neighbor. This reflects the idea that the order in which the connections are drawn is irrelevant.

When considering the ways in which we create connections between nodes, we look to a subset of $\Phi_k[p]$ which is aware of not only the magnitude of a node and the value of $k$, but also the “surroundings” of the node $p$ and the consequences of making the connections indicated by an $\omega$.

Feasible Connection Configurations in $\Gamma_k$.

Define $count(\omega,l)$ for $l \in W$ to be the number of occurrences of symbol $l$ in the word $\omega$. Let $\Gamma_k$ be any $k$-grid and let $p \in \Gamma_k$. The set of all feasible connection configurations is denoted by $\Phi_{\Gamma_k}[p] $ and defined as the set of $\omega \in \Phi_k[p]$ satisfying the following:

* Creating the connections in $\omega$ will not violate the upper bound on connections, $k$.

* Creating the connections in $\omega$ will not create intersecting connections between nodes.

* $count(\omega,l) \leq \min\{k,magn(l_p)\}$ for every $l \in W$.

* The $k$-grid resulting from making the connections in $\omega$ in $\Gamma_k$ does not result in a set of
solved nodes which are disjoint from other nodes.

* No node is incomplete with complete neighbors after making the connections in $\omega$.

Set $\omega^*(p) = \underset{\omega \in \Phi_{\Gamma_k}[p]}{\bigcap} \{ \omega_i \} $ without deleting duplicates.

The conditions required for a word $\omega$ to be in $\Phi_{\Gamma_k}[p]$ follow from the intuition one builds from playing the Numberlink puzzles. The set $\Phi_{\Gamma_k}[p]$ is often a refinement of $\Phi_k[p]$ where now the words in $\Phi_{\Gamma_k}[p]$ obey more of the rules for solving a $k$-grid. Looking at what connections are common to each word in $\Phi_{\Gamma_k}[p]$, gives us a list of connections that are guaranteed to be made in a solution for the $k$-grid.

We show a 2-grid, $\Gamma_2$, on the left with node $p$ in the lower bottom. Applying the word 12 to $p$ we get $\Gamma'_2$ on the right. The node $q$ in $\Gamma'_2$ refers to the center node of magnitude 5.

§.§ Example of computing $\Phi_k[p]$ and $\Phi_{\Gamma_k}[p]$

Let's do an example that illustrates how one constructs the sets $\Phi_k[p]$ and $\Phi_{\Gamma_k}[p]$. Any node of magnitude 2 in a 2-grid will have

$$\Phi_2[p] = \{11, 22, 33, 44, 12, 13, 14, 23, 24, 34\}.$$

Consider the $2$-grid on the left in Figure <ref>.
The node of magnitude 2 in the left bottom corner will have the following as $\Phi_{\Gamma_2}[p]$:

$$\Phi_{\Gamma_2}[p] = \{ 12\}.$$

The words 33, 44, 13, 14, 23, 24, and 34 are not included in the set because of condition 3 in the definition of $\Phi_{\Gamma_2}[p]$ and the words 11 and 22 are not included since they violate condition 4 of the definition of $\Phi_{\Gamma_2}[p]$. We then see that $\omega^*(p)$ for this lower bottom node of magnitude 2 is $\omega^*(p) = 12$. Making these connections gives the $2$-grid on the right of Figure <ref>.

Let's repeat these computations for the center node of magnitude 5, where this time we are working from the 2-grid on the right which includes some connections. We get,

\begin{align*}
\Phi_2[q] = \{11223,11224,11233,11234,11244,11334, & 11344,12233,12234,  12244, \\
& 12334,12344,13344,22334,22344,23344\}
\end{align*}

$$\Phi_{\Gamma'_2}[q] = \{11223, 11224, 11234, 12234\}.$$

In this case, the words that include a 33 will violate condition 3 in the definition of $\Phi_{\Gamma'_2}[q]$ and the words that include a 44 will violate condition 5.
Here we get that $\omega^*(q) = 12$ for the node $q \in \Gamma_2'$. This tells us that any solution to the grid $\Gamma'_2$ will include the connections 12 for node $q$. The reader is encouraged to find a solution to $\Gamma_2$ in Figure <ref>.

§ THE PATH TOWARDS A SOLUTION

In this section we begin with an examination of where to start forming connections in a $k$-grid. We then move to the main result in this work—an algorithm that finds and creates all guaranteed connections between nodes.

§.§ Where do we start?

Let's review the mechanisms we have so far for computing connections between nodes. Based only on the magnitude of a node and $k$ we can generate the list of words in $\Phi_k[n]$. This gives us any and all possible ways to form $n$ connections with 4 or fewer neighbors while respecting the rule that at most $k$ connections are shared between two nodes. We refine this list by looking at the subset $\Phi_{\Gamma_k}[p]$ which is now aware of the neighbors of $p$ and throws away any connection configurations which would immediately result in some violation of the rules for $k$-grids. Finding what is common to all words in $\Phi_{\Gamma_k}[p]$ gives us $\omega^*(p)$.

What is missing thus far is an indication of what node to choose at the start! One way to approach this is to find which nodes have the fewest possible configurations of connections given the node's magnitude, the value of $k$, and the number of neighbors that the node has. For example, a node of magnitude 7 with 4 neighbors in a $2$-grid will have only 4 ways to configure the connections. The large magnitude paired with the small value of $k$ means that there is only 1 connection that is different for each configuration.

\[\includegraphics[scale = .4]{pics/where2start_ex1.png}\]

In Figure <ref> we give a table that shows how many connection configurations exist given a node of magnitude $n$, the value for $k$, and the number of neighbors the node has.
We note some observations based on this table.

* For a fixed $k$ (a row in the table), the number of connection configurations is symmetric as the magnitude $n$ of a node increases. The reader is encouraged to figure out why! Here is a hint: what similarities exist between the connection configurations for the node of magnitude 7 with 4 neighbors and $k=2$ (example above) and a node of magnitude 1 with 4 neighbors and $k=2$?

* The maximum number of connection configurations possible for a fixed $k$ is achieved at magnitude $n = \lfloor\frac{rk}{2}\rfloor$ where $r$ is the number of neighbors  the node has.

* The "tail" of each row in the table for a node with 4 neighbors can be generated by the following formula [8]:

$$ \frac{i(i+1)(i+2)}{6}\hspace{.2in} \text{ for } \hspace{.2in} i =1, 2, 3, \dots$$

* For 3 neighbors, the "tail" of each row is simply the triangular numbers:

$$ \frac{i(i+1)}{2}\hspace{.2in} \text{ for } \hspace{.2in} i =1, 2, 3, \dots.$$

* Where do we start? The table gives us an idea of where to start. The smaller numbers in the table gives us a higher chance of success in creating connections based on $\omega^*$. The table also shows us what nodes are the low-hanging fruit. For example, we see that if a node has $r$ neighbors and has magnitude $magn(p) = rk$, then there is only 1 way to create the needed connections.

Entries in each of the three tables indicate the number of connection configurations that are possible for a node of magnitude $n$ with 4, 3, and 2 neighbors and the indicated $k$. Values in orange indicate when the number of configurations ceases to increase with $k$. Values in blue show what values appear again in the tail of the row below. Values in red are the maximums in each row.

§.§ The $\tau$ algorithm

We now describe an algorithm which creates the guaranteed connections in a $k$-grid. We also show that if this algorithm solves a $k$-grid, then the solution is unique.

To begin we will formalize the idea of creating connections between nodes and how this creates a new $k$-grid.

For any $k$-grid $\Gamma_k$, define the builder function for node $p \in \Gamma_k$ by,

$$u\left(p,\omega^{*}(p)\right) = \Gamma_k^{'}$$

Where $\Gamma_k^{'}$ is $\Gamma_k$ with the following changes:

* $p \in \Gamma_k^{'}$ where $p = \big[x_p, y_p, n_p - ||\omega^{*}||\big]$

* If $q \in \Gamma_k^{'}$ is the $l_p$ neighbor of $p$, then $q = \big[x_q, y_q, n_q - count(l,\omega^{*}) \big]$.

* For each $l \in W$ and pair of nodes $p$ and $l_p$, if $count(l,\omega^{*})>0$, we construct $count(l,\omega^{*})$ connections between $p$ and $l_p$.

At their core, builder functions build connections between a node and its neighbors, reduce the magnitude of the node and reduce the magnitude of the neighboring nodes. We will use these functions to build up a solution (should one exist) for a general $k$-grid.

We would like to apply the builder functions to each node in a $k$-grid again and again in hopes of finding additional guaranteed connections.
We now describe the algorithm which computes and draws the guaranteed connections between nodes. This algorithm will continue until the grid is solved or all guaranteed connections have been exhausted.

$$\textbf{Tau Algorithm}$$

* Choose a node $p$ with $r$ neighbors with $magn(p)=rk$. In this case $k$ connections must be formed with each neighbor. Draw the needed connections for each node $p$ with this property.

* Next, choose a node with one neighbor. Again there is only one destination for those connections! Complete this step for every node $p$ with one neighboring node.

* Now choose a node with one remaining incomplete neighbor. With all other neighbors completed, there is only one place left to draw connecting lines. Draw the needed connections for each node $p$ with one incomplete neighbor.

If there are no remaining nodes with one incomplete neighbor, choose a node $p$ via the following guidelines which are based on the intuition we gather from Figure <ref>. Then compute and draw the connection to create $u(p,\omega^*)$.

* Start with the nodes which have the smallest number of neighbors. Based on Figure <ref>, we know these typically have the fewest possible connection configurations.

* Of these nodes, start with those with the largest and smallest magnitudes (putting off magnitudes that are near $\lfloor \frac{rk}{2} \rfloor$ where $r$ is the number of neighbors that $p$ has).

* Repeat Step 3-4 until either no nodes with magnitude larger than 0 remain or $\omega^*(p) = \emptyset$ for each remaining node with non-zero magnitude.

For any $k$-grid $\Gamma_k$, we denote the resulting grid from the $\tau$ algorithm as $\tau(\Gamma_k)$.

If $\tau(\Gamma_k)$ solves $\Gamma_k$, then $\Gamma_k$ is uniquely solved.

Note that by construction, the word $\omega^*$ for a node $p$ gives a subset of the connections that will exist in any solution to the grid which contains $p$. At each step in the $\tau$ algorithm, we calculate the $\omega^*$ for a node and create the connections indicated by the $\omega^*$. This means that if $\tau(\Gamma_k)$ is a solved $k$-grid, then the set of connections indicated by repeated computations of $\omega^*$ gives the only configuration of connections.  $\blacksquare$

This theorem tells us that if the $\tau$ algorithm solves a $k$-grid, then we can be sure that we have found the only solution possible. One may wonder if there is hope that the above theorem be bidirectional. Unfortunately, it is not true that if a $k$-grid has a unique solution, then the $\tau$ will be able to solve it.  Consider the $k$-grid from Figure <ref>. When $k=2$, this $k$-grid which we call the Pinwheel is uniquely solved.
Every node has $r>1$ neighbors with $magn(p)<rk$, meaning no connections will be drawn from the first three steps. When computing $\omega^*$ for any node we get $\omega^*(p) = \emptyset$. So $\tau$ is unable to make any connections in the Pinwheel!

This $k$-grid which we call the Pinwheel has a unique solution and yet none of it's connections can be found using the $\tau$ algorithm.

A solved $k$-grid whose connections where determined by $\tau$. This $k$-grid is also included in Numberlink and listed as Expert [1].

§ FUTURE WORK

In this work we have created a framework for working with Numberlink puzzles as graphs with certain conditions on the connections between graphs. We have also described an algorithm for computing the guaranteed connections of a node and have shown that in the case that our algorithm solves a $k$-grid, we have found the one unique solution. There are many possible generalizations to this framework and open questions. Below we list a few of these.

* The choice of a square to represent each node is arbitrary and used to emphasize the restriction on the number of neighbors of a node (maximum 4). One could look at graphs whose nodes are polygons (such as an octagon) where connections can be drawn to $m$ neighbors.

* As mentioned in section 3, some nodes have fewer options for how connections can be configured based on $k$ and the number of neighbors the node has. It would be helpful to find a formula to express the number of ways that $n$ connections can be configured for a node based on $k$ and the number of neighbors the node has. This can then tell us what nodes are good starting points for solving the puzzles.
* Having a formula as described above may lead to a more refined algorithm for $\tau$ and may give a way of calculating how many steps will be needed in order to solve a $k$-grid. This can then help find the solution with the least number of steps.

* We know that many $k$-grids are solvable regardless of the value of $k$ and that if a grid is solvable for $k$ then it is also solvable for $k' > k$. Can we find a process to easily compute the smallest value of $k$ such that the grid is solvable?

* In what way can this framework be used for applications in areas like operations research and optimization.

[1]
Puzzledom - numberlink.
<https://play.google.com/store/apps/details?id=com.puzzlegame.puzzledom   hl=en_US   gl=US>.
Accessed: 2021-04-01.

[2]
Daniel J Velleman and Stan Wagon.
Bicycle Or Unicycle?: A Collection of Intriguing Mathematical
Puzzles, volume 36.
MAA Press, 2020.

[3]
Thane Plambeck and Tomas Rokicki.
Barrycades and Septoku: Papers in Honor of Martin Gardner and
Tom Rodgers, volume 100.
American Mathematical Soc., 2020.

[4]
Aaron Adcock, Erik D Demaine, Martin L Demaine, Michael P O'Brien, Felix Reidl,
Fernando Sánchez Villaamil, and Blair D Sullivan.
Zig-zag numberlink is np-complete.
Journal of Information Processing, 23(3):239–245, 2015.

[5]
Suthee Ruangwises and Toshiya Itoh.
Physical zero-knowledge proof for numberlink.
In 10th International Conference on Fun with Algorithms (FUN
2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.

[6]
Ryo Yoshinaka, Toshiki Saitoh, Jun Kawahara, Koji Tsuruma, Hiroaki Iwashita,
and Shin-ichi Minato.
Finding all solutions and instances of numberlink and slitherlink by
Algorithms, 5(2):176–213, 2012.

[7]
Russell Merris.
Graph Theory, volume 1.
Wiley, 2000.

[8]
Neil JA Sloane et al.
The on-line encyclopedia of integer sequences, 2003.


# Synergistic eigenanalysis of covariance and Hessian matrices for enhanced
binary classification

Agus Hartoyo  Jan Argasiński  Aleksandra Trenk  Kinga Przybylska  Anna Błasiak
Alessandro Crimi

###### Abstract

Covariance and Hessian matrices have been analyzed separately in the
literature for classification problems. However, integrating these matrices
has the potential to enhance their combined power in improving classification
performance. We present a novel approach that combines the eigenanalysis of a
covariance matrix evaluated on a training set with a Hessian matrix evaluated
on a deep learning model to achieve optimal class separability in binary
classification tasks. Our approach is substantiated by formal proofs that
establish its capability to maximize between-class mean distance and minimize
within-class variances. By projecting data into the combined space of the most
relevant eigendirections from both matrices, we achieve optimal class
separability as per the linear discriminant analysis (LDA) criteria. Empirical
validation across neural and health datasets consistently supports our
theoretical framework and demonstrates that our method outperforms established
methods. Our method stands out by addressing both LDA criteria, unlike PCA and
the Hessian method, which predominantly emphasize one criterion each. This
comprehensive approach captures intricate patterns and relationships,
enhancing classification performance. Furthermore, through the utilization of
both LDA criteria, our method outperforms LDA itself by leveraging higher-
dimensional feature spaces, in accordance with Cover’s theorem, which favors
linear separability in higher dimensions. Our method also surpasses kernel-
based methods and manifold learning techniques in performance. Additionally,
our approach sheds light on complex DNN decision-making, rendering them
comprehensible within a 2D space.

###### keywords:

covariance matrix , Hessian matrix , eigenanalysis , binary classification ,
class separability

††journal: ar$\chi$iv

[first]organization=Sano - Centre for Personalised Computational
Medicine,city=Krakow, country=Poland [third]organization=Jagiellonian
University,city=Krakow, country=Poland [second]organization=School of
Computing, Telkom University,city=Bandung, country=Indonesia

##  1 Introduction

Binary classification is a fundamental task in machine learning, where the
goal is to assign data points to one of two classes. The accuracy and
effectiveness of binary classifiers depend on their ability to separate the
two classes accurately. However, achieving optimal class separability can be
challenging, especially when dealing with complex and high-dimensional data.

Traditional approaches often rely on analyzing either the covariance matrix
[1, 2, 3, 4, 5, 6, 7] or the Hessian matrix [8, 9, 10, 11, 12, 13, 14]
separately to optimize machine learning models. In the field of Evolution
Strategies (ESs), a recent work [15] explored the relationship between the
covariance matrix and the landscape Hessian, highlighting the statistical
learning capabilities of ESs using isotropic Gaussian mutations and rank-based
selection. While this study provides valuable insights into ESs’ learning
behavior, it does not investigate the practical integration of the two
matrices.

When analyzing the covariance matrix, the focus is on capturing the inherent
patterns of variability within the data, allowing for a compact representation
that highlights relationships between different dimensions. On the other hand,
the analysis of the Hessian matrix aims to find the direction in which the
classes are best separated, by maximizing the curvature along the
discriminative directions. However, these separate analyses fail to fully
leverage the synergistic effects that can arise from integrating the
information contained in both matrices.

To tackle this challenge, we present a novel approach that combines the
eigenanalysis of the covariance matrix evaluated on a training set with the
Hessian matrix evaluated on a deep learning model. By integrating these
matrices, our method aims to optimize class separability by simultaneously
maximizing between-class mean distance and minimizing within-class variances,
which are the fundamental criteria of linear discriminant analysis (LDA) [16,
17]. By utilizing these criteria, we leverage the foundational principles that
have been established and proven over decades of LDA’s application in various
fields like medicine [18, 19, 20, 21, 22], agriculture [23, 24, 25, 26, 27],
and biometrics [28, 29, 30, 31, 32, 33, 34]. LDA has been demonstrating the
reliability and usefulness of the two criteria as indicators of discriminative
power. Thus, by building upon this well-established foundation, our approach
will inherit the strengths and reliability that have been demonstrated by LDA.
This integrated approach holds the promise of enhancing classification
performance beyond conventional methods that treat the two matrices
separately.

##  2 Methodology

###  2.1 Approach: combining eigenanalysis of covariance and Hessian matrices

In this subsection, we describe our novel approach for combining the
eigenanalysis of the covariance matrix and the Hessian matrix to achieve
optimal class separability in binary classification tasks.

  1. 1.

Covariance matrix eigenanalysis: A covariance matrix
$\text{Cov}(\boldsymbol{\theta})$ is a $D\times D$ matrix, where $D$
represents the dimensionality of the predictor attributes. Each element of the
covariance matrix reflects the covariance between the two predictor attributes
$\theta_{1}$ and $\theta_{2}$:

$\text{{Cov}}(\theta_{1},\theta_{2})=\frac{1}{n-1}\sum_{i=1}^{n}(\theta_{1i}-\bar{\theta}_{1})(\theta_{2i}-\bar{\theta}_{2})$

where $\theta_{1i}$ and $\theta_{2i}$ are the corresponding observations of
these predictor attributes for the $i$-th data instance, $\bar{\theta}_{1}$
and $\bar{\theta}_{2}$ are the sample means of the predictor attributes
$\theta_{1}$ and $\theta_{2}$, respectively, and $n$ is the number of data
instances or observations.

Performing eigenanalysis on the covariance matrix
$\text{Cov}(\boldsymbol{\theta})$, we obtain eigenvalues $\lambda_{i}$ and
corresponding eigenvectors $\mathbf{v}_{i}$. The leading eigenvector
$\mathbf{v}_{1}$ associated with the largest eigenvalue $\lambda_{1}$, as
captured in the eigen-equation
$\text{Cov}(\boldsymbol{\theta})\cdot\mathbf{v}_{1}=\lambda_{1}\cdot\mathbf{v}_{1}$,
represents the principal direction with the highest variance.

  2. 2.

Hessian matrix eigenanalysis: Next, we compute a Hessian matrix evaluated on a
deep learning model trained on the same training set. We employ a deep neural
network (DNN) architecture consisting of four fully connected layers, each
followed by a Rectified Linear Unit (ReLU) activation. The final layer employs
a sigmoid activation function to yield a probability value within the range of
0 to 1. During training, we utilize the binary cross-entropy loss function.
The binary cross-entropy loss is given by

$\text{BCELoss}=-\sum_{i=1}^{n}\log p_{\theta}(c_{i}\mid x_{i}),$

where $c_{i}$ is the actual class for data instance $x_{i}$. The probability
$p_{\theta}(c_{i}\mid x_{i})$ is the predicted probability of instance $x_{i}$
belonging to class $c_{i}$, as estimated by the DNN model parameterized by
$\theta$.

The Hessian of the loss function is then:

$\displaystyle H_{\theta}$ $\displaystyle=\nabla_{\theta}^{2}\text{BCELoss}$
$\displaystyle=\nabla_{\theta}^{2}\left[-\sum_{i=1}^{n}\log
p_{\theta}(c_{i}\mid x_{i})\right]$
$\displaystyle=-\sum_{i=1}^{n}\left[\nabla_{\theta}^{2}\log
p_{\theta}(c_{i}\mid x_{i})\right].$

Performing eigenanalysis on the Hessian matrix $H_{\boldsymbol{\theta}}$, we
obtain eigenvalues $\lambda_{i}^{\prime}$ and corresponding eigenvectors
$\mathbf{v}_{i}^{\prime}$. The leading eigenvector $\mathbf{v}_{1}^{\prime}$
associated with the largest eigenvalue $\lambda_{1}^{\prime}$, as captured in
the eigen-equation
$H_{\boldsymbol{\theta}}\cdot\mathbf{v}_{1}^{\prime}=\lambda_{1}^{\prime}\cdot\mathbf{v}_{1}^{\prime}$,
represents the direction corresponding to the sharpest curvature.

  3. 3.

Integration of matrices and projection of data:

To combine the power of covariance and Hessian matrices, we project the data
into the combined space of the most relevant eigendirections. Let $\mathbf{U}$
be the matrix containing the leading eigenvectors from both matrices:

$\mathbf{U}=[\mathbf{v}_{1},\mathbf{v}_{1}^{\prime}]$ (1)

The 2D projection of the data, denoted as $\mathbf{X}_{\text{proj}}$, is
obtained by:

$\mathbf{X}_{\text{proj}}=\mathbf{X}\cdot\mathbf{U}$ (2)

Here, $\mathbf{X}$ is the original data matrix, and $\mathbf{X}_{\text{proj}}$
represents the final output of the proposed method—a 2D projection capturing
both statistical spread and discriminative regions of the data.

Unlike LDA, which aims to optimize both criteria simultaneously along a single
direction for binary classification, constrained by the limitation that the
number of linear discriminants is at most $c-1$ where $c$ is the number of
class labels [35], our approach provides more flexibility and control. By
working on two separate directions, we specifically focus on minimizing the
within-class variances in one direction while maximizing the between-class
mean distance in the other direction.

###  2.2 Formal foundation: maximizing the squared between-class mean distance
and minimizing the within-class variance

In this subsection, we establish two theorems along with their respective
proof sketches that form the theoretical basis for our approach. Full proofs
for both theorems are available in Appendix A.

####  2.2.1 Theorem 1: Maximizing Covariance for Maximizing Squared Between-
Class Mean Distance.

Consider two sets of 1D data points representing two classes, denoted as
$C_{1}$ and $C_{2}$, each consisting of $n$ samples. The data in $C_{1}$ and
$C_{2}$ follow the same underlying distribution centered around their
respective means, $m_{1}$ and $m_{2}$. Here, $C_{2}$ can be understood as a
reflection of $C_{1}$ with the axis of reflection positioned at the overall
mean, denoted as $m$. Furthermore, the variances of $C_{1}$ and $C_{2}$, are
equal, denoted as $s_{1}^{2}=s_{2}^{2}=s_{w}^{2}$, and the combined data from
$C_{1}$ and $C_{2}$ has a variance of $s^{2}$. The between-class mean
distance, denoted as $d$, represents the separation between the means of
$C_{1}$ and $C_{2}$. We establish the following relationship:

$s^{2}=\frac{d^{2}}{4(1-\lambda)}$ (3)

where $\lambda=\frac{s_{w}^{2}}{s^{2}}$ signifies a constant between 0 and 1
that reflects the distribution of the original data when $C_{1}$ and $C_{2}$
are considered as projected representations.

##### Proof sketch:

  1. 1.

Consider two classes $C_{1}$ and $C_{2}$ with identical underlying
distributions.

  2. 2.

Simplify the expression for the combined data variance $s^{2}$ to
$s_{w}^{2}+\frac{1}{4}d^{2}$.

  3. 3.

Apply the Variance Ratio Preservation Theorem for Projection onto a Vector to
relate $s_{w}^{2}$ and $s^{2}$ through $\lambda$.

  4. 4.

Simplify the expression to obtain $s^{2}=\frac{d^{2}}{4(1-\lambda)}$.

  5. 5.

Conclude that maximizing $s^{2}$ maximizes $d^{2}$, fulfilling the theorem’s
objective.

####  2.2.2 Theorem 2: Maximizing Hessian for Minimizing Within-Class Variance

Let $\theta$ be a parameter of the model, and $\mathrm{H}_{\theta}$ denote the
Hessian of the binary cross-entropy loss function with respect to $\theta$. We
define the within-class variance as the variance of a posterior distribution
$p_{\theta}(\theta\mid c_{i})$, which represents the distribution of the
parameter $\theta$ given a class $c_{i}$. We denote the variance of this
posterior distribution as $\sigma_{post}^{2}$. We establish the following
relationship:

$\mathrm{H}_{\theta}\approx\frac{1}{\sigma_{\text{post}}^{2}}.$ (4)

##### Proof sketch:

  1. 1.

Consider the parameter $\theta$ and the Hessian $\mathrm{H}_{\theta}$ of the
binary cross-entropy loss function. Define the within-class variance
$\sigma_{\text{post}}^{2}$ as the variance of the posterior distribution
$p_{\theta}(\theta\mid c_{i})$.

  2. 2.

Approximate the Hessian $\mathrm{H}_{\theta}$ as the Fisher information using
the expectation of the squared gradient of the log-likelihood [36].

  3. 3.

Assume a known normal likelihood distribution for
$p_{\theta}(c_{i}\mid\theta)$ with mean $\mu$ and standard deviation $\sigma$.
Compute the Fisher information as $\frac{1}{\sigma^{2}}$.

  4. 4.

Considering a uniform prior distribution $p(\theta)$ within the plausible
range of $\theta$ and the evidence $p(c_{i})$ as a known constant, apply
Bayes’ formula to derive that $\sigma_{\text{post}}^{2}=\sigma^{2}$.

  5. 5.

Derive $\mathrm{H}_{\theta}\approx\frac{1}{\sigma_{\text{post}}^{2}}$.

  6. 6.

Conclude that maximizing $\mathrm{H}_{\theta}$ minimizes
$\sigma_{\text{post}}^{2}$, achieving the theorem’s objective.

Theorem 1 and Theorem 2 respectively suggest that maximizing the variance of
projected data and the Hessian effectively maximize squared between-class mean
distance and minimize within-class variances. These theorems provide the
theoretical foundation for the eigenanalysis of the covariance and Hessian
matrices, crucial steps in our proposed method for improving class
separability based on the two LDA criteria.

##  3 Complexity analysis

In assessing the computational demands of our proposed method, we employ big O
notation to describe the worst-case time complexity. The overall time
complexity of our method encompasses several tasks, including covariance and
Hessian matrix computations ($\mathcal{O}(N\cdot F^{2})$), eigenanalysis of
these matrices ($\mathcal{O}(F^{3})$), selection of eigenvectors corresponding
to leading eigenvalues ($\mathcal{O}(1)$), and data projection into the
combined space ($\mathcal{O}(N\cdot F)$). Here, $N$ signifies the number of
data points in the training set, while $F$ represents the number of features
per data point. It is essential to note that this analysis considers our
method’s computational demands under the assumption that a pre-existing DNN is
in place.

Comparatively, the time complexity of LDA is primarily influenced by the
calculation of both within-class and between-class scatter matrices
($\mathcal{O}(N\cdot F^{2})$) and subsequent eigenanalysis
($\mathcal{O}(F^{3})$). In both methods, the dominant complexity factor is
either the matrix computations ($\mathcal{O}(N\cdot F^{2})$) when there are
more examples than features or the eigenanalysis step ($\mathcal{O}(F^{3})$)
otherwise.

This analysis reveals that our proposed method exhibits a comparable
computational profile to the method under improvement. Therefore, our approach
offers enhanced class separability and interpretability without significantly
increasing computational demands, making it a practical choice for real-world
applications.

##  4 Experiment

###  4.1 Assessed dimensionality reduction techniques

In our evaluation, we compare nine distinct dimensionality reduction and data
projection techniques, each offering unique insights: principal component
analysis (PCA), kernel PCA (KPCA) [37], Hessian, uniform manifold
approximation and projection (UMAP) [38], locally linear embedding (LLE) [39],
linear optimal low-rank projection (LOL) [40], linear discriminant analysis
(LDA), kernel discriminant analysis (KDA) [41], and the proposed approach. PCA
involves projection onto the primary two covariance eigenvectors, i.e., it
applies KPCA with a linear kernel. Similarly, LDA employs KDA with a linear
kernel. The Hessian method projects data onto the leading two Hessian
eigenvectors. Notably, for both KPCA and KDA, the best kernels (linear,
polynomial, RBF, sigmoid, or cosine similarity) and their associated
parameters (kernel coefficient or degree) are determined using grid search
tailored to each dataset. The proposed method combines the most relevant
eigenvectors derived from both the covariance and Hessian matrices, as
specified by Eqs(1, 2).

###  4.2 Class separability assessment via linear SVMs

The nine dimensionality reduction methods we assess are designed to transform
high-dimensional data into more manageable 1D or 2D spaces while
simultaneously enhancing or preserving class separability. Linear SVMs enable
us to create decision boundaries that are readily visualized in these 1D or 2D
spaces. Thus, we utilize SVMs with a linear kernel to evaluate and visualize
the extent of class separability achieved through the projections.

In addition to the visualization of decision boundaries, we employ a
comprehensive set of evaluation metrics to quantify the performance of the
dimensionality reduction methods. These metrics include the F1 score, which
measures the balance between precision and recall, the Area Under the Receiver
Operating Characteristic Curve (AUC ROC), which assesses the model’s ability
to distinguish between positive and negative classes, and Cohen’s kappa, a
statistic that gauges the agreement between the predicted and actual class
labels.

###  4.3 Datasets

We will conduct our assessment of the nine distinct methods on the following
datasets:

Widely recognized benchmark datasets. We evaluate our approach using three
widely recognized benchmark datasets for binary classification: the Wisconsin
breast cancer dataset [42], the heart disease dataset [43], and the Pima
Indians diabetes dataset [44]. Prior to applying various dimensionality
reduction methods, we enact standard data preprocessing techniques on the
original datasets, including handling of missing data, one-hot encoding for
categorical variables, and data normalization.

Neural spike train dataset. A spike train is a sequence of action potentials
(spikes) emitted by a neuron over time. The neural spike train dataset used in
this research consists of recordings from rat’s neurons during drug
application from a multi-electrode array [45, 46]. The data, comprising 221
records, represents the final dataset after all preprocessing steps suggested
by [47]. Each record contains 15 time-series features extracted using the
‘tsfresh‘ package [48] and 1 class attribute indicating whether the neuron is
non-responsive (0) or responsive (1). Exploratory data analysis revealed an
imbalance in the dataset, with 190 non-responsive neurons (86%) and 31
responsive neurons (14%).

To provide an unbiased assessment of our method’s performance, we conduct
extensive experiments on the datasets with 10-fold cross-validation with the
exception of the neural spike train dataset. Due to the limited number of
positive cases in the neural spike train dataset, we performed 5-fold cross-
validation to ensure a reasonable sample size for validation.

###  4.4 Results

Figure 1: Projection of the Wisconsin breast cancer data into different
combined spaces of the covariance and Hessian eigenvectors. (a) Nine selected
projection plots, each representing data projected onto a distinct space
created by combining the first three covariance and first three Hessian
eigenvectors. (b) Heatmap showing the squared between-class mean distance for
projections onto varying combinations of covariance and Hessian eigenvectors.
The heatmap demonstrates that the values remain constant vertically across
different Hessian eigenvectors, while exhibiting a noticeable descending order
horizontally, aligning with the descending order of the variance. These
results essentially concretize our formal premise, empirically validating the
linear relationship described in Eq(3) between the variance and the squared
between-class mean distance. (c) Heatmap showing the sum of within-class
variances for projections onto different combinations of covariance and
Hessian eigenvectors. The values remain constant horizontally across different
covariance eigenvectors but exhibit a clear ascending order vertically,
aligning with the descending order of the Hessian. These empirical results
validate the negative correlation between the Hessian and the within-class
variance described in Eq(4) within the framework of our theoretical
foundation. (d) Heatmap displaying the LDA ratio, representing the ratio
between the squared between-class mean distances presented in (b) and the
corresponding within-class variances shown in (c). The highest LDA ratio is
observed for the combination of the first Hessian eigenvector with the first
covariance eigenvector. Notably, a general descending pattern is observed both
horizontally and vertically across different combinations, indicating that
both covariance eigenanalysis (represented along the horizontal direction) and
Hessian eigenanalysis (represented along the vertical direction) equally
contribute to the class separability (represented by the LDA ratio).

The results depicted in Figure 1, which pertain to the WBCD dataset, along
with the corresponding findings presented in Appendix B for the other three
datasets, collectively reaffirm the consistency and robustness of our proposed
approach across diverse datasets. The heatmaps shown in Figure 1(b) and
Appendix B consistently demonstrate the reduction in squared between-class
mean distance as more covariance eigenvectors are incorporated, aligning with
our established theoretical framework (Eq(3)). Furthermore, Figure 1(c) and
its counterparts in the appendix reveal the ascending order of within-class
variances, in sync with the descending order of the Hessian, supporting our
theoretical foundation (Eq(4)). Furthermore, Figure 1(d) and its counterparts
illustrate the LDA ratio, emphasizing the equal contributions of covariance
and Hessian eigenanalyses to class separability. Importantly, the results also
imply that the highest LDA ratio is observed for the combination of the first
Hessian eigenvector with the first covariance eigenvector. This observation
underscores the significance of projecting data onto the combined space of
these primary eigenvectors (as outlined in Eqs (1, 2)), forming the core of
our proposed method. These consistent empirical results across multiple
datasets not only validate our theoretical premises but also endorse the
effectiveness of our proposed method in optimizing class separability.

The evaluation results in Figure 2 compare the performance of different data
projection methods, including PCA, KPCA, Hessian, UMAP, LLE, LOL, LDA, KDA,
and the proposed method, using 5- or 10-fold cross-validation. It is essential
to note that we introduced nonlinearity to the comparative analysis by
utilizing distinct kernels determined through grid search in KPCA and KDA for
each dataset. Notably, our proposed method consistently outperforms all
others, securing the highest scores across all datasets and evaluation
metrics. This consistent superiority of our approach implies its potential as
a valuable tool for improving classification performance in various domains,
further highlighting its promise as a robust and effective method for
dimensionality reduction and data projection.

Figure 2: Performance comparison of data projection methods using cross-
validation on four distinct datasets: (a) WBCD, (b) heart disease, (c) neural
spike train, and (d) Pima Indians diabetes datasets. This figure presents the
average F1 score, ROC AUC, and Cohen’s Kappa values obtained through 5- or
10-fold cross-validation for nine data projection techniques: PCA, KPCA
(cosine similarity for WBCD, polynomial kernel with degree=3 for Heart, linear
kernel for Neural, cosine similarity for Pima), Hessian, UMAP, LLE, LOL, LDA,
KDA (cosine similarity for WBCD, RBF kernel with coefficient=0.01 for Heart,
sigmoid kernel with coefficient=2 for Neural, cosine similarity for Pima), and
the proposed method. Notably, the proposed method consistently outperforms all
other techniques, achieving the highest scores across all evaluation metrics.

##  5 Discussion

Our work provides a compelling theoretical insight and a powerful, practical
method, demonstrating the strength of simplicity in achieving remarkable
results. Firstly, we provide a profound theoretical insight, revealing a
subtle yet powerful relationship between covariance and Hessian matrices. Our
formal proof seamlessly links covariance eigenanalysis with the first LDA
criterion while Hessian eigenanalysis with the second one. This unification
under LDA criteria offers a fresh and intuitive perspective on their
interplay. Secondly, capitalizing on this theoretical elegance and simplicity,
we introduce a novel method that consistently outperforms established
techniques across diverse datasets. The unexpected efficacy of our method,
rooted in the straightforward relationship between covariance, Hessian, and
LDA, showcases the effectiveness of simplicity in addressing complex
challenges.

The proposed method outperforms PCA and the Hessian method by comprehensively
addressing both LDA criteria—maximizing between-class mean distance and
minimizing within-class variances. Unlike PCA, which predominantly focuses on
the former and lacks the guidance of class labels, our supervised approach
considers both aspects. Despite the computational efficiency associated with
unsupervised dimension reduction methods [49], our approach demonstrates the
added value of incorporating class labels. This key insight also underlies our
outperformance of KPCA, which, despite operating on non-linearities, remains
essentially unsupervised in nature. While the Hessian method concentrates on
minimizing within-class variances, our method optimally combines the strengths
of PCA and the Hessian, effectively identifying feature space directions that
enhance both between-class separation and within-class compactness.

While the proposed method employs the LDA criteria, it surpasses LDA itself in
all cases. Figure 3 visually demonstrates the advantages of the proposed
method over LDA. LDA is limited to a one-dimensional projection for binary
classification problems [35], where it seeks to identify a single direction
that simultaneously satisfies the two criteria. Conversely, the proposed
method splits the task of meeting the criteria into two directions. The
utilization of higher dimensionality in the proposed method increases the
likelihood of discovering class separability, aligning with Cover’s theorem
[50]. KDA, operating on a non-linear mode, outperforms LDA in all cases, yet
it remains fundamentally confined to one dimension, restricting its
effectiveness in capturing intricate class-specific patterns compared to our
proposed method.

Figure 3: Projection and classification results of UMAP, LDA, and the proposed
method with SVM on the Wisconsin breast cancer test data. This figure
showcases the test data projected by three distinct data projection methods,
each separated by its corresponding SVM hyperplane trained using the WBCD
training data. This figure highlights the effectiveness of linear SVM in
facilitating class separability visualization and interpretability of the
model.

Our method outperforms LOL by leveraging non-linear modes of operation,
providing a distinct advantage in capturing complex patterns beyond the linear
capabilities of LOL. UMAP exhibits good class separability with the widest
margin as shown in Figure 3. However, it is important to note that UMAP is not
inherently a classification technique, and thus, it fails to generalize well
to new data. Similar limitations apply to LLE, which, although effective in
revealing local data structures, lacks the inherent capability for
classification and generalization.

Figure 3 also underscores the simplicity and interpretability of linear SVMs
as basic linear classifiers in dealing with low-dimensional data. The figure
shows clear visualizations of the SVM’s decision boundaries and separation
achieved through the various dimensionality reduction methods. Notably, our
proposed method combined with linear SVMs offers valuable insights into the
decision-making process within the underlying DNN. Leveraging the transparency
and interpretability provided by SVMs to address the inherent opacity of DNNs,
our method bridges the gap between complex DNNs and comprehensible, highly
accurate decision processes within a 2D space.

The proposed method shows promising performance, but it has certain
limitations. The applicability of our approach is based on the premise that a
DNN model is already in place. Our method operates on top of the underlying
DNN. Any shortcomings or biases in the DNN’s performance will naturally
reflect in the results obtained from our approach. The limitation of our
method lies in its dependence on the quality and accuracy of the underlying
DNN.

An important trajectory for future work involves investigating the extension
of our methodology to accommodate various loss functions beyond binary cross-
entropy. The mathematical derivation in our current work relies on the elegant
relationship between binary cross-entropy loss and within-class variances.
Exploring the adaptability of our method to different loss functions will
contribute to a more comprehensive understanding of the method’s versatility,
but requires careful scrutiny to establish analogous connections.
Simultaneously, we recognize the need to extend our methodology from binary to
multiclass classification. The binary classification focus in this work stems
from foundational aspects guiding our formal proof, which is designed around
binary assumptions to facilitate a streamlined and elegant derivation process.
In particular, the use of binary cross-entropy as the loss function and the
utilization of a linear SVM for evaluation inherently adhere to binary
classification. Moving forward, careful exploration is needed to adapt our
approach to multiclass scenarios to ensure its applicability and effectiveness
across a broader range of classification tasks.

##  6 Conclusion

Our paper presents a multifaceted contribution to the realms of binary
classification and dimensionality reduction. We offer a rigorous formal proof,
showcasing how our novel method, grounded in the eigenanalysis of covariance
and Hessian matrices, systematically enhances class separability, drawing
inspiration from LDA criteria. This theoretical foundation finds strong
empirical support through comprehensive experiments across an array of
datasets. In these experiments, our approach consistently outperforms
established techniques, demonstrating its robustness and applicability.
Moreover, our approach shares a similar complexity profile with the
traditional method, ensuring its practical utility in real-world scenarios.
Additionally, through the integration of our method with linear SVMs, we
improve the explainability of the intricate decision-making processes inherent
to DNNs, addressing their typical lack of transparency and facilitating
enhanced model interpretability.

##  7 Reproducibility statement

The detailed proofs for the theoretical foundations, emphasizing the
maximization of between-class mean distance and minimization of within-class
variance, are provided in Appendix A. The complete source code and datasets
for our experiments are accessible through the following links:

  1. 1.

WBCD dataset experiment: Link to Colab notebook

  2. 2.

Heart disease dataset experiment: Link to Colab notebook

  3. 3.

Neural spike train dataset experiment: Link to Colab notebook

  4. 4.

Pima Indians diabetes dataset experiment: Link to Colab notebook

These notebooks contain the complete source code along with the datasets
accessed from the same Google Drive account, facilitating easy reproduction
and comprehension of our results.

## References

  * [1] T.-H. Nagai, The covariance matrix of green’s functions and its application to machine learning, ArXiv abs/2004.06481 (2020).
  * [2] H. Q. Minh, V. Murino, Covariances in computer vision and machine learning, in: Covariances in Computer Vision and Machine Learning, 2017.
  * [3] G. Serra, C. Grana, M. Manfredi, R. Cucchiara, Covariance of covariance features for image classification, Proceedings of International Conference on Multimedia Retrieval (2014).
  * [4] K. Lenc, A. Vedaldi, Learning covariant feature detectors, in: ECCV Workshops, 2016.
  * [5] P. D. Hoff, X. Niu, A covariance regression model, Statistica Sinica 22 (2011) 729–753.
  * [6] B.-C. Kuo, D. A. Landgrebe, A covariance estimator for small sample size classification problems and its application to feature extraction, IEEE Trans. Geosci. Remote. Sens. 40 (2002) 814–819.
  * [7] C. Lam, High‐dimensional covariance matrix estimation, Wiley Interdisciplinary Reviews: Computational Statistics 12 (2019).
  * [8] A. Dawid, P. Huembeli, M. Tomza, M. Lewenstein, A. Dauphin, Hessian-based toolbox for reliable and interpretable machine learning in physics, Machine Learning: Science and Technology 3 (2021).
  * [9] S. Fu, W. Liu, D. Tao, Y. Zhou, L. Nie, Hesgcn: Hessian graph convolutional networks for semi-supervised classification, Inf. Sci. 514 (2020) 484–498.
  * [10] Z. Yao, A. Gholami, K. Keutzer, M. W. Mahoney, Pyhessian: Neural networks through the lens of the hessian, 2020 IEEE International Conference on Big Data (Big Data) (2019) 581–590.
  * [11] G. Krishnasamy, R. Paramesran, Hessian semi-supervised extreme learning machine, Neurocomputing 207 (2016) 560–567.
  * [12] S. Wiesler, J. Li, J. Xue, Investigations on hessian-free optimization for cross-entropy training of deep neural networks, in: Interspeech, 2013.
  * [13] R. H. Byrd, G. M. Chin, W. Neveitt, J. Nocedal, On the use of stochastic hessian information in optimization methods for machine learning, SIAM J. Optim. 21 (2011) 977–995.
  * [14] J. Martens, Deep learning via hessian-free optimization, in: International Conference on Machine Learning, 2010.
  * [15] O. M. Shir, A. Yehudayoff, On the covariance-hessian relation in evolution strategies, Theoretical Computer Science 801 (2020) 157–174.
  * [16] R. A. Fisher, The use of multiple measurements in taxonomic problems, Annals of eugenics 7 (2) (1936) 179–188.
  * [17] P. Xanthopoulos, P. M. Pardalos, T. B. Trafalis, P. Xanthopoulos, P. M. Pardalos, T. B. Trafalis, Linear discriminant analysis, Robust data mining (2013) 27–33.
  * [18] A. Sharma, S. Imoto, S. Miyano, A between-class overlapping filter-based method for transcriptome data analysis, Journal of bioinformatics and computational biology 10 (05) (2012) 1250010.
  * [19] A. Sharma, K. K. Paliwal, Cancer classification by gradient lda technique using microarray gene expression data, Data & Knowledge Engineering 66 (2) (2008) 338–347.
  * [20] B. Moghaddam, Y. Weiss, S. Avidan, Generalized spectral bounds for sparse lda, in: Proceedings of the 23rd international conference on Machine learning, 2006, pp. 641–648.
  * [21] S. Dudoit, J. Fridlyand, T. P. Speed, Comparison of discrimination methods for the classification of tumors using gene expression data, Journal of the American statistical association 97 (457) (2002) 77–87.
  * [22] H.-P. Chan, D. Wei, M. A. Helvie, B. Sahiner, D. D. Adler, M. M. Goodsitt, N. Petrick, Computer-aided classification of mammographic masses and normal tissue: linear discriminant analysis in texture feature space, Physics in Medicine & Biology 40 (5) (1995) 857.
  * [23] A. Tharwat, T. Gaber, A. E. Hassanien, One-dimensional vs. two-dimensional based features: Plant identification approach, Journal of Applied logic 24 (2017) 15–31.
  * [24] T. Gaber, A. Tharwat, V. Snasel, A. E. Hassanien, Plant identification: Two dimensional-based vs. one dimensional-based feature extraction methods, in: 10th international conference on soft computing models in industrial and environmental applications, Springer, 2015, pp. 375–385.
  * [25] S. Rezzi, D. E. Axelson, K. Héberger, F. Reniero, C. Mariani, C. Guillou, Classification of olive oils using high throughput flow 1h nmr fingerprinting with principal component analysis, linear discriminant analysis and probabilistic neural networks, Analytica Chimica Acta 552 (1-2) (2005) 13–24.
  * [26] K. Héberger, E. Csomós, L. Simon-Sarkadi, Principal component and linear discriminant analyses of free amino acids and biogenic amines in hungarian wines, Journal of agricultural and food chemistry 51 (27) (2003) 8055–8060.
  * [27] L. Chen, N. C. Carpita, W.-D. Reiter, R. H. Wilson, C. Jeffries, M. C. McCann, A rapid method to screen for cell-wall mutants using discriminant analysis of fourier transform infrared spectra, The Plant Journal 16 (3) (1998) 385–392.
  * [28] K. K. Paliwal, A. Sharma, Improved pseudoinverse linear discriminant analysis method for dimensionality reduction, International Journal of Pattern Recognition and Artificial Intelligence 26 (01) (2012) 1250002.
  * [29] L. Yuan, Z.-c. Mu, Ear recognition based on 2d images, in: 2007 First IEEE International Conference on Biometrics: Theory, Applications, and Systems, IEEE, 2007, pp. 1–5.
  * [30] C. H. Park, H. Park, Fingerprint classification using fast fourier transform and nonlinear discriminant analysis, Pattern Recognition 38 (4) (2005) 495–503.
  * [31] X. Wang, X. Tang, Random sampling lda for face recognition, in: Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004., Vol. 2, IEEE, 2004, pp. II–II.
  * [32] H. Yu, J. Yang, A direct lda algorithm for high-dimensional data—with application to face recognition, Pattern recognition 34 (10) (2001) 2067–2070.
  * [33] L.-F. Chen, H.-Y. M. Liao, M.-T. Ko, J.-C. Lin, G.-J. Yu, A new lda-based face recognition system which can solve the small sample size problem, Pattern recognition 33 (10) (2000) 1713–1726.
  * [34] R. Haeb-Umbach, H. Ney, Linear discriminant analysis for improved large vocabulary continuous speech recognition., in: icassp, Vol. 92, Citeseer, 1992, pp. 13–16.
  * [35] J. Ye, R. Janardan, Q. Li, Two-dimensional linear discriminant analysis, Advances in neural information processing systems 17 (2004).
  * [36] E. Barshan, M.-E. Brunet, G. K. Dziugaite, Relatif: Identifying explanatory training samples via relative influence, in: International Conference on Artificial Intelligence and Statistics, PMLR, 2020, pp. 1899–1909.
  * [37] B. Schölkopf, A. Smola, K.-R. Müller, Kernel principal component analysis, in: International conference on artificial neural networks, Springer, 1997, pp. 583–588.
  * [38] L. McInnes, J. Healy, J. Melville, Umap: Uniform manifold approximation and projection for dimension reduction, arXiv preprint arXiv:1802.03426 (2018).
  * [39] S. T. Roweis, L. K. Saul, Nonlinear dimensionality reduction by locally linear embedding, science 290 (5500) (2000) 2323–2326.
  * [40] J. T. Vogelstein, E. W. Bridgeford, M. Tang, D. Zheng, C. Douville, R. Burns, M. Maggioni, Supervised dimensionality reduction for big data, Nature communications 12 (1) (2021) 2872.
  * [41] S. Mika, G. Ratsch, J. Weston, B. Scholkopf, K.-R. Mullers, Fisher discriminant analysis with kernels, in: Neural networks for signal processing IX: Proceedings of the 1999 IEEE signal processing society workshop (cat. no. 98th8468), Ieee, 1999, pp. 41–48.
  * [42] W. N. Street, W. H. Wolberg, O. L. Mangasarian, Nuclear feature extraction for breast tumor diagnosis, in: Biomedical image processing and biomedical visualization, Vol. 1905, SPIE, 1993, pp. 861–870.
  * [43] R. Detrano, A. Janosi, W. Steinbrunn, M. Pfisterer, J.-J. Schmid, S. Sandhu, K. H. Guppy, S. Lee, V. Froelicher, International application of a new probability algorithm for the diagnosis of coronary artery disease, The American journal of cardiology 64 (5) (1989) 304–310.
  * [44] J. W. Smith, J. E. Everhart, W. Dickson, W. C. Knowler, R. S. Johannes, Using the adap learning algorithm to forecast the onset of diabetes mellitus, in: Proceedings of the annual symposium on computer application in medical care, American Medical Informatics Association, 1988, p. 261.
  * [45] D. Tsai, E. John, T. Chari, R. Yuste, K. Shepard, High-channel-count, high-density microelectrode array for closed-loop investigation of neuronal networks, in: 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2015, pp. 7510–7513.
  * [46] M. O. Heuschkel, M. Fejtl, M. Raggenbass, D. Bertrand, P. Renaud, A three-dimensional multi-electrode array for multi-site stimulation and recording in acute brain slices, Journal of neuroscience methods 114 (2) (2002) 135–148.
  * [47] I. Lazarevich, I. Prokin, B. Gutkin, V. Kazantsev, Spikebench: An open benchmark for spike train time-series classification, PLOS Computational Biology 19 (1) (2023) e1010792.
  * [48] M. Christ, N. Braun, J. Neuffer, A. W. Kempa-Liehr, Time series feature extraction on basis of scalable hypothesis tests (tsfresh–a python package), Neurocomputing 307 (2018) 72–77.
  * [49] C. Shen, M. Sun, M. Tang, C. E. Priebe, Generalized canonical correlation analysis for classification, Journal of Multivariate Analysis 130 (2014) 310–322.
  * [50] T. M. Cover, Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition, IEEE transactions on electronic computers (3) (1965) 326–334.

##  Appendix A Full proofs: maximizing the squared between-class mean distance
and minimizing the within-class variance

###  A.1 Proof of Theorem 1 - Maximizing Covariance for Maximizing Squared
Between-Class Mean Distance

To prove that maximizing the variance will maximize the squared between-class
mean distance, we start by considering two sets of 1D data points representing
two classes, denoted as $C_{1}$ and $C_{2}$, each consisting of $n$ samples.
The data in $C_{1}$ and $C_{2}$ follow the same underlying distribution
centered around their respective means, $m_{1}$ and $m_{2}$. In other words,
$C_{2}$ can be understood as a reflection of $C_{1}$ with the axis of
reflection positioned at the overall mean, denoted as $m$, effectively giving
rise to a shifted configuration. Specifically, the indices for $C_{1}$ range
from 1 to $n$, while the indices for $C_{2}$ range from $n+1$ to $2n$.

For any pair of data points $x_{i}$ and $x_{j}$, where $1\leq i\leq n$ and
$n+1\leq j\leq 2n$, the difference between $x_{i}$ and $m_{1}$ is equal to the
difference between $x_{j}$ and $m_{2}$, i.e., $x_{i}-m_{1}=x_{j}-m_{2}$. This
relationship ensures that the shifted copies of $C_{1}$ and $C_{2}$ maintain
the same relative distances from their respective means, preserving the
identical distribution in both sets. Furthermore, the variances of $C_{1}$ and
$C_{2}$, are equal, denoted as $s_{1}^{2}=s_{2}^{2}=s_{w}^{2}$, and the
combined data from $C_{1}$ and $C_{2}$ has a variance of $s^{2}$.

The between-class mean distance, denoted as $d$, represents the separation
between the means of $C_{1}$ and $C_{2}$. We can express the means as
$m_{1}=m-\frac{d}{2}$ and $m_{2}=m+\frac{d}{2}$, where $m$ is effectively
located in the middle, equidistant from both $m_{1}$ and $m_{2}$.

The variance of the combined data is given by:

$s^{2}=\frac{1}{2n-1}\left(\sum_{i=1}^{n}(x_{i}-m)^{2}+\sum_{i=n+1}^{2n}(x_{i}-m)^{2}\right)$

To simplify the expression, we substitute $m=m_{1}+\frac{1}{2}d$ and
$m=m_{2}-\frac{1}{2}d$. This yields:

$\displaystyle s^{2}$
$\displaystyle=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(x_{i}-\left(m_{1}+\frac{1}{2}d\right)\right)^{2}+\sum_{i=n+1}^{2n}\left(x_{i}-\left(m_{2}-\frac{1}{2}d\right)\right)^{2}\right)$
$\displaystyle=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)-\frac{1}{2}d\right)^{2}+\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)+\frac{1}{2}d\right)^{2}\right)$
$\displaystyle\begin{split}&=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)^{2}-\left(x_{i}-m_{1}\right)d+\frac{1}{4}d^{2}\right)\right)\\\
&+\left.\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)^{2}+\left(x_{i}-m_{2}\right)d+\frac{1}{4}d^{2}\right)\right)\\\
&=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)^{2}+\frac{1}{4}d^{2}\right)-\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)d\right)\right)\\\
&+\left.\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)^{2}+\frac{1}{4}d^{2}\right)\right)+\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)d\right)\\\
&=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)^{2}+\frac{1}{4}d^{2}\right)+\left.\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)^{2}+\frac{1}{4}d^{2}\right)\right)\right)\\\
&+\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)d\right)-\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)d\right)\\\
&=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)^{2}+\frac{1}{4}d^{2}\right)+\left.\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)^{2}+\frac{1}{4}d^{2}\right)\right)\right)\\\
&+\left(\sum_{i=n+1}^{2n}\left(x_{i}-m_{2}\right)-\sum_{i=1}^{n}\left(x_{i}-m_{1}\right)\right)d\\\
\end{split}$

The property of the identical distribution in both sets suggests
$\sum_{i=n+1}^{2n}\left(x_{i}-m_{2}\right)-\sum_{i=1}^{n}\left(x_{i}-m_{1}\right)=0$.
So, we get:

$\displaystyle s^{2}$
$\displaystyle=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(\left(x_{i}-m_{1}\right)^{2}+\frac{1}{4}d^{2}\right)+\sum_{i=n+1}^{2n}\left(\left(x_{i}-m_{2}\right)^{2}+\frac{1}{4}d^{2}\right)\right)$
$\displaystyle=\frac{1}{2n-1}\left(\sum_{i=1}^{n}\left(x_{i}-m_{1}\right)^{2}+\sum_{i=n+1}^{2n}\left(x_{i}-m_{2}\right)^{2}+\frac{1}{2}n\cdot
d^{2}\right)$
$\displaystyle\approx\frac{1}{2\left(n-1\right)}\left(\sum_{i=1}^{n}\left(x_{i}-m_{1}\right)^{2}+\sum_{i=n+1}^{2n}\left(x_{i}-m_{2}\right)^{2}\right)+\frac{1}{4}d^{2}$
$\displaystyle=\frac{1}{2}\left(\frac{\sum_{i=1}^{n}\left(x_{i}-m_{1}\right)^{2}}{\left(n-1\right)}+\frac{\sum_{i=n+1}^{2n}\left(x_{i}-m_{2}\right)^{2}}{\left(n-1\right)}\right)+\frac{1}{4}d^{2}$
$\displaystyle=\frac{1}{2}\left(s_{w}^{2}+s_{w}^{2}\right)+\frac{1}{4}d^{2}$
$\displaystyle=s_{w}^{2}+\frac{1}{4}d^{2}.$

Now, considering the data points representing the projected data onto an
(Eigen)vector, we can utilize the Variance Ratio Preservation Theorem for
Projection onto a Vector, which establishes the relationship between $s^{2}$
and $s_{w}^{2}$ as follows:

$\displaystyle s_{w}^{2}=\lambda\cdot s^{2}$

where $\lambda$ is a constant between 0 and 1, determined by the distribution
of the original data being projected.

Substituting this equation into the previous expression, we have:

$\displaystyle s^{2}$ $\displaystyle=\lambda\cdot s^{2}+\frac{1}{4}d^{2}.$

Let’s rearrange the equation by moving $2\lambda\cdot s^{2}$ to the left side:

$\displaystyle s^{2}-\lambda\cdot s^{2}$ $\displaystyle=\frac{1}{4}d^{2}.$

Combining like terms:

$\displaystyle(1-\lambda)\cdot s^{2}$ $\displaystyle=\frac{1}{4}d^{2}.$

To solve for $s^{2}$, divide both sides by $(1-\lambda)$:

$\displaystyle s^{2}$
$\displaystyle=\frac{\frac{1}{4}d^{2}}{1-\lambda}=\frac{r^{2}}{1-\lambda}$

where $r=\frac{1}{2}d=m-m_{1}=m_{2}-m$.

We observe that the sign of $s^{2}$ and $d^{2}$ will be the same since the
denominator $1-\lambda$ is always positive (as $0<\lambda<1$). Therefore,
$s^{2}$ is linearly proportional to $d^{2}$.

Hence, maximizing the variance ($s^{2}$) will maximize the squared between-
class mean distance ($d^{2}$) as desired.

###  A.2 Proof of Theorem 2 - Maximizing Hessian for Minimizing Within-Class
Variance

We aim to prove that maximizing the Hessian will minimize the within-class
variance. Let $\theta$ denote a parameter of the classifier.

We define the within-class variance as the variance of a posterior
distribution $p_{\theta}(\theta\mid c_{i})$, which represents the distribution
of the parameter $\theta$ given a class $c_{i}$. We denote the variance of
this posterior distribution as $\sigma_{post}^{2}$.

Recall that our Hessian is given by:

$\displaystyle\mathrm{H}_{\theta}$
$\displaystyle=-\left[\nabla_{\theta}^{2}\log
p_{\theta}(c_{i}\mid\theta)\right]$

However, according to [36], we can approximate the Hessian using Fisher
information:

$\mathrm{H}_{\theta}\approx\mathbb{E}_{p_{\theta}}\left[\left(\nabla_{\theta}\log
p_{\theta}(c_{i}\mid\theta)\right)^{2}\right]$ (5)

Assuming that a known normal distribution underlies the likelihood
$p_{\theta}(c_{i}\mid\theta)$, i.e.,

$p_{\theta}(c_{i}\mid\theta)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(\theta-\mu)^{2}}{2\sigma^{2}}\right)$

where $\mu$ and $\sigma$ are known constants, we can compute the Fisher
information as follows:

$\displaystyle\mathbb{E}_{p_{\theta}}\left[\left(\nabla_{\theta}\log
p_{\theta}(c_{i}\mid\theta)\right)^{2}\right]$
$\displaystyle=\mathbb{E}_{p_{\theta}}\left[\left(\nabla_{\theta}\log\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(\theta-\mu)^{2}}{2\sigma^{2}}\right)\right)^{2}\right]$
$\displaystyle=\mathbb{E}_{p_{\theta}}\left[\left(\frac{\theta-\mu}{\sigma^{2}}\right)^{2}\right]$
$\displaystyle=\frac{1}{\sigma^{4}}\mathbb{E}_{p_{\theta}}\left[(\theta-\mu)^{2}\right]$
$\displaystyle=\frac{1}{\sigma^{4}}\sigma^{2}$
$\displaystyle=\frac{1}{\sigma^{2}}$

From Eq(5), we derive the Hessian as follows:

$\mathrm{H}_{\theta}\approx\frac{1}{\sigma^{2}}$ (6)

Now, let us assume that the evidence $p(c_{i})$ is a known constant $\rho$. We
also assume that the prior distribution $p(\theta)$ follows a uniform
distribution within the plausible range of $\theta$, which is bounded by a
known minimum value $\theta_{\min}$ and maximum value $\theta_{\max}$.
Formally:

$p_{\theta}(\theta)=\begin{cases}\dfrac{1}{\theta_{\max}-\theta_{\min}},&\text{if}\
\theta_{\min}\leq\theta\leq\theta_{\max}\\\ 0,&\text{otherwise}\end{cases}$

Based on Bayes’ formula and the given assumptions, the posterior distribution
within the plausible range of $\theta$ is:

$\displaystyle p_{\theta}(\theta\mid c_{i})$
$\displaystyle=\frac{p_{\theta}(c_{i}\mid\theta)\cdot
p_{\theta}(\theta)}{p(c_{i})}$
$\displaystyle=\frac{1}{\rho(\theta_{\max}-\theta_{\min})}\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(\theta-\mu)^{2}}{2\sigma^{2}}\right)$

This implies that the posterior distribution $p_{\theta}(\theta\mid c_{i})$ is
a normal distribution with mean $\mu$ and variance $\sigma^{2}$:

$p_{\theta}(\theta\mid c_{i})\sim\mathcal{N}(\mu,\,\sigma^{2})$

Therefore, the within-class variance $\sigma_{\text{post}}^{2}$ of the
posterior distribution is equal to $\sigma^{2}$.

Combining the previous result with the Hessian calculation in Eq(6), we can
conclude that:

$\displaystyle\mathrm{H}_{\theta}$
$\displaystyle\approx\frac{1}{\sigma_{\text{post}}^{2}}$

Hence, maximizing the Hessian matrix ($\mathrm{H}_{\theta}$) will minimize the
within-class variances ($\sigma_{\text{post}}^{2}$) as desired.

###  A.3 Supporting theorems

####  A.3.1 Distance Preservation Theorem for Projection onto a Vector

In a 2-dimensional space with a first axis and a second axis, consider a set
of points located solely on the first axis, denoted by $(x_{i},0)$, where
$x_{i}$ represents the position of a point along the first axis. Let
$\mathbf{v}=(v_{1},v_{2})$ be a unit vector, with $v_{1}$ and $v_{2}$ as its
components. The projection operation maps each point $(x_{i},0)$ to its
corresponding scalar projection $y_{i}$ onto the vector $\mathbf{v}$. The
Distance Preservation Theorem states that the projection onto $\mathbf{v}$
preserves the ratio of distances between the points on the first axis.

##### Proof:

The scalar projection $y_{i}$ of a point $(x_{i},0)$ onto the vector
$\mathbf{v}$ is given by:

$y_{i}=\frac{(x_{i},0)\cdot\mathbf{v}}{|\mathbf{v}|}=\frac{(x_{i},0)\cdot\mathbf{v}}{1}=x_{i}\cdot
v_{1}+0\cdot v_{2}=x_{i}\cdot v_{1}$

This implies that the projection of each point $(x_{i},0)$ from the first axis
onto the vector $\mathbf{v}$ is obtained by multiplying the coordinate $x_{i}$
by the first component $v_{1}$ of $\mathbf{v}$.

Now, considering the distances between two points $(x_{i},0)$ and $(x_{j},0)$
on the first axis and their corresponding projections $y_{i}$ and $y_{j}$ onto
the vector $\mathbf{v}$, we define the distance between $(x_{i},0)$ and
$(x_{j},0)$ as:

$d_{1}=|x_{i}-x_{j}|$

Similarly, the distance between $y_{i}$ and $y_{j}$ is defined as:

$d_{2}=|y_{i}-y_{j}|$

Substituting the expressions for $y_{i}$ and $y_{j}$ derived earlier, we have:

$d_{2}=|x_{i}\cdot v_{1}-x_{j}\cdot v_{1}|=|v_{1}|\cdot|x_{i}-x_{j}|$

It can be observed that the ratio of the distances is constant:

$\frac{d_{2}}{d_{1}}=\frac{|v_{1}|\cdot|x_{i}-x_{j}|}{|x_{i}-x_{j}|}=|v_{1}|$

This shows that the ratio of distances between the points on the first axis is
preserved in the projection onto $\mathbf{v}$.

Therefore, the Distance Preservation Theorem for Projection onto a Vector
concludes that the projection operation onto the vector
$\mathbf{v}=(v_{1},v_{2})$ in a 2-dimensional space preserves the ratio of
distances between the points on the first axis.

####  A.3.2 Variance Preservation Theorem for Projection onto a Vector

Consider an arbitrary subset of points located on the first axis, denoted by
$X=\\{x_{1},x_{2},\ldots,x_{n}\\}$, with a mean $\mu_{X}$ and variance
$\sigma_{X}^{2}$. These points are projected onto the vector
$\mathbf{v}=(v_{1},v_{2})$ using the projection operation defined earlier.

The projected subset of points on $\mathbf{v}$ is denoted by
$Y=\\{y_{1},y_{2},\ldots,y_{n}\\}$, where each $y_{i}$ represents the
projection of $x_{i}$ onto $\mathbf{v}$. The variance of $Y$, denoted as
$\sigma_{Y}^{2}$, is a measure of the spread of the projected points around
its mean $\mu_{Y}$.

The Variance Preservation Theorem states that the ratio of variances between
the projected subset $Y$ and the original subset $X$ is equal to the square of
the first component $v_{1}$ of the projection vector $\mathbf{v}$.

Mathematically, this can be expressed as:

$\frac{\sigma_{Y}^{2}}{\sigma_{X}^{2}}=v_{1}^{2}$

##### Proof:

The variance of $Y$ can be computed as:

$\sigma_{Y}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(y_{i}-\mu_{Y})^{2}$

Similarly, the variance of $X$ can be computed as:

$\sigma_{X}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\mu_{X})^{2}$

According to the Distance Preservation Theorem for Projection onto a Vector,
the squared distances between the points on the first axis and their
projections can be related as:

$(y_{i}-\mu_{Y})^{2}=|y_{i}-\mu_{Y}|^{2}=(|v_{1}|\cdot|x_{i}-\mu_{X}|)^{2}=v_{1}^{2}\cdot(x_{i}-\mu_{X})^{2}$

Substituting this expression into the variance of $Y$, we have:

$\sigma_{Y}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}v_{1}^{2}\cdot(x_{i}-\mu_{X})^{2}=v_{1}^{2}\cdot\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\mu_{X})^{2}=v_{1}^{2}\cdot\sigma_{X}^{2}$

Therefore, we have shown that:

$\sigma_{Y}^{2}=v_{1}^{2}\cdot\sigma_{X}^{2}$

Hence, the ratio of variances between the projected subset $Y$ and the
original subset $X$ is given by $v_{1}^{2}$, as stated in the Variance
Preservation Theorem for Projection onto a Vector.

####  A.3.3 Variance Ratio Preservation Theorem for Projection onto a Vector

Consider an arbitrary subset of points located on the first axis, denoted by
$X=\\{x_{1},x_{2},\ldots,x_{n}\\}$, with variances $\sigma_{X_{1}}^{2}$ and
$\sigma_{X_{2}}^{2}$ for subsets $X_{1}$ and $X_{2}$, respectively. These
points are projected onto the vector $\mathbf{v}=(v_{1},v_{2})$ using the
projection operation defined earlier.

The projected subsets of points on $\mathbf{v}$ are denoted by
$Y_{1}=\\{y_{1},y_{2},\ldots,y_{n}\\}$ and
$Y_{2}=\\{z_{1},z_{2},\ldots,z_{n}\\}$, where each $y_{i}$ and $z_{i}$
represents the projection of $x_{i}$ onto $\mathbf{v}$. The variances of
$Y_{1}$ and $Y_{2}$, denoted as $\sigma_{Y_{1}}^{2}$ and $\sigma_{Y_{2}}^{2}$
respectively, measure the spread of the projected points around their
respective means.

The Variance Ratio Preservation Theorem states that the ratio of variances
between the projected subsets $Y_{2}$ and $Y_{1}$ is equal to the ratio of
variances between the original subsets $X_{2}$ and $X_{1}$.

Mathematically, this can be expressed as:

$\frac{\sigma_{Y_{2}}^{2}}{\sigma_{Y_{1}}^{2}}=\frac{\sigma_{X_{2}}^{2}}{\sigma_{X_{1}}^{2}}$

##### Proof:

By the Variance Preservation Theorem for Projection onto a Vector, we know
that the variance of the projected subset $Y$ can be expressed as:

$\sigma_{Y}^{2}=v_{1}^{2}\cdot\sigma_{X}^{2}$

Applying this theorem to $Y_{1}$ and $X_{1}$, we have:

$\sigma_{Y_{1}}^{2}=v_{1}^{2}\cdot\sigma_{X_{1}}^{2}$

Similarly, applying the theorem to $Y_{2}$ and $X_{2}$, we have:

$\sigma_{Y_{2}}^{2}=v_{1}^{2}\cdot\sigma_{X_{2}}^{2}$

Now, we can take the ratio of these variances:

$\frac{\sigma_{Y_{2}}^{2}}{\sigma_{Y_{1}}^{2}}=\frac{v_{1}^{2}\cdot\sigma_{X_{2}}^{2}}{v_{1}^{2}\cdot\sigma_{X_{1}}^{2}}=\frac{\sigma_{X_{2}}^{2}}{\sigma_{X_{1}}^{2}}$

Therefore, we have shown that the ratio of variances between the projected
subsets $Y_{2}$ and $Y_{1}$ is equal to the ratio of variances between the
original subsets $X_{2}$ and $X_{1}$, as required.

##  Appendix B Additional dataset results

Figure 4: Projection of the heart disease data into different combined spaces
of the covariance and Hessian eigenvectors. Figure 5: Projection of the neural
spike train data into different combined spaces of the covariance and Hessian
eigenvectors. Figure 6: Projection of the Pima Indians diabetes data into
different combined spaces of the covariance and Hessian eigenvectors. (a) Nine
selected projection plots, each representing data projected onto a distinct
space created by combining the first three covariance and first three Hessian
eigenvectors. (b) Heatmap showing the squared between-class mean distance for
projections onto varying combinations of covariance and Hessian eigenvectors.
(c) Heatmap showing the sum of within-class variances for projections onto
different combinations of covariance and Hessian eigenvectors. (d) Heatmap
displaying the LDA ratio, representing the ratio between the squared between-
class mean distances presented in (b) and the corresponding within-class
variances shown in (c).