# Iterative Methods at Lower Precision

Yizhou Chen, Xiaoyun Gong, Xiang Ji
Guest Editor: James<EMAIL_ADDRESS>
###### Abstract

Since numbers in the computer are represented with a fixed number of bits,
loss of accuracy during calculation is unavoidable. At high precision where
more bits (e.g. 64) are allocated to each number, round-off errors are
typically small. On the other hand, calculating at lower precision, such as
half (16 bits), has the advantage of being much faster. This research focuses
on experimenting with arithmetic at different precision levels for large-scale
inverse problems, which are represented by linear systems with ill-conditioned
matrices. We modified the Conjugate Gradient Method for Least Squares (CGLS)
and the Chebyshev Semi-Iterative Method (CS) with Tikhonov regularization to
do arithmetic at lower precision using the MATLAB chop function, and we ran
experiments on applications from image processing and compared their
performance at different precision levels. We concluded that CGLS is a more
stable algorithm, but overflows easily due to the computation of inner
products, while CS is less likely to overflow but it has more erratic
convergence behavior. When the noise level is high, CS outperforms CGLS by
being able to run more iterations before overflow occurs; when the noise level
is close to zero, CS appears to be more susceptible to accumulation of round-
off errors.

##  1 Introduction

Most computer processors today use double-precision binary floating point
arithmetic, which represents floating-point numbers with 64 bits. The
convention follows from the IEEE-754 standard established in 1985, specifying
the number of bits for the sign, exponent, and mantissa (fraction) for various
floating-point formats including binary16 (half precision), binary32 (single
precision), and binary64 (double precision). The difference between these
formats is illustrated in the diagram below [4]:

Taking a quarter of memory of their traditional 64-bit counterparts, computing
with half precision numbers has attracted interest from researchers and
manufacturers because of the potential for significantly reduced computational
time. The format accelerates deep learning training, allows larger models, and
improves gaming latency, providing players with better experiences [12].
However, aside from these benefits, half precision also brings one obvious
disadvantage: increased inaccuracies in number representation. One goal of
this research project is to investigate the impact of transferring from double
precision operations to low precision arithmetic when solving large-scale
inverse problems.

Due to their wide variety of applications, we specifically utilized and
modified codes for iterative methods in solving inverse problems to evaluate
the performance of low precision algorithms. Inverse problems arise when using
outside measurements to acquire information about internal or hidden data [6].
We operate on known outputs with some errors to compute the true input value.
For example, X-ray computed tomography is a contactless imaging method for
reconstructing target objects, most commonly, pathologies in the human body.
Another example is image deblurring problems, which occur when the true
picture is to be reconstructed from its blurry (and sometimes noisy) version
[2]. These cases can be abstracted as linear systems
${\bm{b}}=A{\bm{x}}+{\bm{e}}$, where $A$ is a large-scale, typically ill-
conditioned matrix and ${\bm{b}}$ is a vector output blended with noise,
${\bm{e}}$.

A naive solution to this problem can be obtained by directly calculating a
solution to $A{\bm{x}}={\bm{b}}$. However, this naive solution is often
corrupted by noise due to the ill-conditioning of matrix $A$. Regularization
methods are often needed to balance signal and noise. One approach for
regularization uses the singular value decomposition of $A$. However, when the
matrix $A$ is large, as is often the case for many real-life applications
(including our test cases), such direct methods may be difficult to implement,
since decomposing the matrix can be very computationally costly. When the
matrix is sparse, meaning it contains a lot of zero entries, computing matrix-
vector products can be done very efficiently. Therefore, for such problems,
implementing iterative methods, which only require matrix-vector products and
vector operations, is a more efficient way to solve inverse problems than
direct methods, with the former requiring $O(n^{3})$ work for an $n\times n$
matrix $A$, while the latter involves significantly less than $O(n^{2})$ work
per iteration [13].

##  2 Simulating Low Precision

###  2.1 Chop

The MATLAB function chop written by Higham and Pranesh provides users a way to
simulate low precision numbers and arithmetic. The function can simulate
precisions, such as fp16, bfloat16, as well as custom formats where the user
is able to choose the number of bits in the significand and the maximum value
of the exponent. The function doesn’t create the new data type; instead, it
keeps numbers in double or single precision and makes sure that they have the
right scale as in lower precision. Therefore, once the rounding meets the bit
limit of the corresponding target precision, the remaining bits are set to
zero. Computer operations are also carried out in double or single precision,
and after that, the result is rounded to low precision [10].

Fully simulating computations in low precision requires users to call chop
after each arithmetic operation. For example, for the operation

$a=x+y\times z,$

the appropriate use of chop is:

$a=chop(x+chop(y\times z)).$

This is burdensome yet unavoidable, so our modified version of iterative
methods using chop, which only simulates low precision arithmetic, usually
takes a long time to run, despite the fact that low precision arithmetic would
run significantly faster with the corresponding hardware. However, for vector
or matrix operations, there are ways to reduce the number of calls to chop,
improving the efficiency of our codes. For example, for the vector inner
product between ${\bm{x}}$ and ${\bm{y}}$, instead of using the line:

    sum = 0;
    for i = 1:vector_length
        sum = chop(sum + chop(x(i) * y(i)));
    end

we use an element-wise operation at the beginning:

    sum = 0;
    z = chop(x.*y)
    for i = 1:vector_length
        sum = chop(sum + z(i));
    end

such that we successfully reduce $2n$ calls of chop to only $n+1$ calls,
accelerating our running process.

###  2.2 Blocking

Since we are using floating-point arithmetic, inaccuracy is inevitable in
computations. However, blocking can be used to reduce the error bound. The
method breaks a large number of operations into several smaller pieces,
computes them independently, and sums them up. Consider the inner product
between two vectors:

${\bm{x}}=[x_{1},x_{2},x_{3},x_{4},x_{5},x_{6}],\;\;{\bm{y}}=[y_{1},y_{2},y_{3},y_{4},y_{5},y_{6}]$

Instead of calculating it directly, we could break it into:

${\bm{x}}_{1}=[x_{1},x_{2},x_{3}],\;\;{\bm{x}}_{2}=[x_{4},x_{5},x_{6}]$
${\bm{y}}_{1}=[y_{1},y_{2},y_{3}],\;\;{\bm{y}}_{2}=[y_{4},y_{5},y_{6}]$

Then we calculate ${\bm{x}}_{1}^{T}{\bm{y}}_{1}$ and
${\bm{x}}_{2}^{T}{\bm{y}}_{2}$, and add the sum. The result can have less
errors than the direct calculation. The error bound for inner products
${\bm{x}}^{T}{\bm{y}}$ is [8]:

$|{\bm{x}}^{T}{\bm{y}}-fl({\bm{x}}^{T}{\bm{y}})|\leq\gamma_{n}|{\bm{x}}|^{T}|{\bm{y}}|$

where $\displaystyle\gamma_{n}=\frac{nu}{1-nu}$ and $u$ is the unit round-off
of floating point arithmetic (e.g. $9.77*10^{-4}$ for half precision,
$1.19*10^{-7}$ for single precision, and $2.22*10^{-16}$ for double
precision). Nonetheless, with blocking, the error bound is reduced to [8]:

$|s_{n}-\hat{s}_{n}|\leq\gamma_{(\log_{2}n)+1}|{\bm{x}}|^{T}|{\bm{y}}|.$

$s_{n}$ is the summation of the results for each block using the true value of
${\bm{x}}$ and ${\bm{y}}$, whereas $\hat{s}_{n}$ is the summation of each
block’s result using the floating number representation of ${\bm{x}}$ and
${\bm{y}}$. The detailed proof is in [8, Chapter 3]. We plotted the error
graph for double, single, and half precision against block size when computing
the inner product of two random vectors in Figure 1. We did inner products for
each precision $20$ times and calculated the average error, which is computed
as the difference between the result from our chopped version of inner
products and MATLAB’s default double precision computation.

Figure 1: Relationship between errors and block size.

The graph on the left-hand-side includes all three precision levels, and we
can see that the error for half precision is the largest, since it has the
least bits. Focusing on half precision only, the graph on the right shows that
errors decrease sharply as blocking is introduced. However, the errors
increase as the block size increases, indicating that an appropriate block
size needs to be carefully chosen. The reason behind this increase is that as
the block size grows larger, it has the similar effect as doing no blocking;
the whole vector or matrix is put into the first block, not divided into
smaller sections. In our modified codes, we chose the block size to be $256$,
appropriate for the $4096\times 4096$ matrices we used in the test cases.

##  3 Conjugate Gradient Method for Least Squares

###  3.1 Method Overview

The Conjugate Gradient (CG) Method was introduced by Hestenes and Stiefel [7].
It is an iterative method for solving the linear system $A{\bm{x}}={\bm{b}}$,
where $A$ is a symmetric and positive definite matrix. The CG method can be
viewed as an optimization problem of minimizing a convex quadratic function:

$\phi({\bm{x}})=\frac{1}{2}{\bm{x}}^{T}A{\bm{x}}-{\bm{x}}^{T}{\bm{b}}.$

The gradient is zero at minimum point, meaning
$\nabla\phi({\bm{x}})=A{\bm{x}}-{\bm{b}}=0$, which is exactly the linear
system we are trying to solve.

The CG method is a Krylov subspace method, which means its approximated
solutions lie in Krylov subspaces. In each iteration, ${\bm{x}}$ is allowed to
explore subspaces of increasing dimensions. An interesting property of CG is
that during each iteration, ${\bm{x}}$ is updated to the point within the
subspaces where the $A$-norm of the error is minimized [13]. If we assume zero
round-off errors, CG is guaranteed to converge within a limited number of
steps. Specifically, if $A$ is an $n\times n$ matrix, the method will find the
solution in no greater than $n$ steps. Besides, the updated ${\bm{x}}$ is
closer to the true solution than the previous ${\bm{x}}$ in each iteration.

The CGLS algorithm is the least squares version of the CG method, applied to
the normal equations $A^{T}A{\bm{x}}=A^{T}{\bm{b}}$.

One potential problem with this method for low precision is that the
calculation of inner products can easily result in overflow, as illustrated in
the experiment section below.

The following algorithm describes CGLS [1].

Algorithm 1 Conjugate Gradient Method Least Squares

Let ${\bm{x}}^{(0)}$ be an initial approximation, set ${\bf
r}^{(0)}={\bm{b}}-A{\bm{x}}^{(0)},{\bf p}^{(0)}={\bf s}^{(0)}=A^{T}{\bf
r}^{(0)},\psi_{0}=||{\bf s}^{(0)}||_{2}^{2}$

while $\psi_{k}>tol$ do

${\bf q}^{(k)}=A{\bf p}^{(k)}$,

$\alpha_{k}=\psi_{k}/||{\bf q}^{(k)}||_{2}^{2}$,

${\bm{x}}^{(k+1)}={\bm{x}}^{(k)}+\alpha_{k}{\bf p}^{(k)}$,

${\bf r}^{(k+1)}={\bf r}^{(k)}+\alpha_{k}{\bf q}^{(k)}$,

${\bf s}^{(k+1)}=A^{T}{\bf r}^{(k+1)}$,

$\psi_{k+1}=||{\bf s}^{(k+1)}||_{2}^{2}$,

$\beta_{k}=\psi_{k+1}/\psi_{k}$,

${\bf p}^{(k+1)}={\bf s}^{(k+1)}+\beta_{k}{\bf p}^{(k)}$.

end while

###  3.2 Experiment

With the chop function, we modified the original CGLS method in the IRtools
package, which offers various iterative methods for large-scale, ill-posed
inverse problems and a set of test problems for these iterative methods [3].
Then we used two large-scale, ill-posed inverse problems from the same
package: image deblurring and tomography reconstruction. The first is to
reconstruct an approximation of the true image from the observed blurry
version, whereas the second one is to reconstruct an image from measured
projections, which can be obtained, for example, from X-ray beams. We
investigated these two test problems in different sizes, floating-point
precision levels, and noise levels. In each of these examples, the problem is
modeled as ${\bm{b}}=A{\bm{x}}+{\bm{e}}$, where $A$ and ${\bm{b}}$ are given.

####  3.2.1 Image Deblurring

Here we generate an image deblurring test problem using the PRblur function in
the IRtools package. The true image is a picture of the Hubble space
telescope, and the observed data is corrupted by Gaussian blurring function;
for details, see [3]. We first added no noise to ${\bm{b}}$ (i.e.
${\bm{e}}=0$) to see how our modified CGLS method works. In Figures 4, 4, and
4, we displayed the computed approximation of ${\bm{x}}$ at the final
iteration.

Figure 2: Double precision, size 64, zero noise.

Figure 3: Single precision, size 64, zero noise.

Figure 4: Half precision, size 64, zero noise.

The single-precision result is very similar to that computed in double-
precision, but the half-precision result is more blurry, and the background
contains more artifacts. We also plotted the error graphs for all three
formats in Figure 5.

Figure 5: The error norm of a size 64 problem with mild blurring of different
precisions.

All three error norms overlap from the beginning until around $20^{th}$
iteration, where the half-precision errors begins to deviate from those in
single and double precision. The difference is due to the round-off errors of
half precision, which add up and take over. Moreover, the error norms of half
precision terminates at $28^{th}$ iteration, because overflow of inner
products causes NaNs (Not a Number) to be computed during the iteration.

In real life, the noise-free image, ${\bm{b}}$, is not often obtainable but
most likely it is blended with noise. Therefore, it is imperative to make sure
that the chopped CGLS method works effectively with noisy data as well. We
added noise at different levels to ${\bm{b}}$ and displayed the computed
approximations of ${\bm{x}}$ in the last iteration in half precision Figures
8, 8, and 8.

Figure 6: Half precision, size 64, 0.1% noise.

Figure 7: Half precision, size 64, 1% noise.

Figure 8: Half precision, size 64, 10% noise.

However, the last iteration doesn’t necessarily mean it’s the iteration with
the best performance, so the results generated using ${\bm{x}}$ from the best
iteration (i.e. the iteration with smallest error norms) are also shown in
Figures 11, 11, and 11.

Figure 9: Half precision, problem size 64 with mild blurring and 0.1% noise at
best iteration.

Figure 10: Half precision, problem size 64 with mild blurring and 1% noise at
best iteration.

Figure 11: Half precision, problem size 64 with mild blurring and 10% noise at
best iteration.

For results in the last iteration, while images obtained when using half
precision are more distorted than their single counterparts, they follow the
same trend in terms of the impact of noise. $0.1\%$ noise behaves very
similarly to the cases with no noise; however, $1\%$ noise begins to dominate,
and the background has substantially more artifacts, while the object is still
identifiable. Eventually, the computed result is horribly corrupted by
artifacts in the 10% noise case; the picture contains very little meaningful
information.

Figure 12: Error norm for problem size 64 with mild blurring and 0.1% noise.

Figure 13: Error norm for problem size 64 with mild blurring and 1% noise.

Figure 14: Error norm for problem size 64 with mild blurring and 10% noise.

An interesting phenomenon emerges if we examine the error norm for each
iteration in Figures 14, 14, and 14. While the norm decreases as the iteration
continues for the $0.1\%$ noise cases, those with $1\%$ and $10\%$ noise
present a different trend: the error norm starts to increase again after
certain iterations. The reason for this reversal is that the noise in
${\bm{b}}$ is inverted during each iteration along with the true ${\bm{b}}$;
the noisy data eventually outweigh the actual solution, and the error norm
increases. This phenomenon is known as semi-convergence in the inverse
problems literature [6].

So far all experimental results are acquired using our chopped version of the
CGLS code without regularization. Section $4$ discusses another iterative
method called Chebyshev semi-iterative method, which has the advantage of not
requiring any inner products. However, the Chebyshev method does not have the
same convergence behavior as CGLS, and it is therefore necessary to
incorporate another regularization technique known as Tikhonov regularization.
In order to make comparisons between these two methods, we later added
Tikhonov regularization in our modified CGLS. The results and graphs will be
displayed in Section $5$.

####  3.2.2 Tomography Reconstruction

Tomography reconstruction is another type of inverse problem that produces
images from X-ray projection data of various angles. We used it as a test
problem since the IRtools software package includes a simulation in the PRtomo
function. The modified CGLS method works well for double and single precision
and generates nice reconstructions as shown in Figures 17 and 17, but issues
arise for half precision because entries of the solution ${\bm{x}}$ are all
NaNs. Since half precision formats allocate only $5$ bits for the exponent,
overflow easily occurs when the inner product is calculated. NaNs are the
results from dividing infinity by infinity in the CGLS method.

One of the solutions to the overflow is to rescale the matrices. After scaling
both $A$ and ${\bm{b}}$ by $0.01$, we obtain Figure 17. Although the image is
still not as clear as those from double and single precision, the inner
products did not overflow and a meaningful picture was obtained.

Figure 15: Double precision, size 64, zero noise.

Figure 16: Single precision, size 64, zero noise.

Figure 17: Half precision, size 64, zero noise (after rescaling).

The noise level affects the quality of the reconstruction as well. As Figures
20, 20, 20 illustrate, while random noise does not completely destroy the
reconstruction and take over the original information as is the case for image
deblurring, the artifacts caused by the noise make it difficult to see small
objects in the image.

Figure 18: Half precision, size 64, zero noise.

Figure 19: Half precision, size 64, 1% noise.

Figure 20: Half precision, size 64, 10% noise (after rescaling).

As was the case with the deblurring problem, the error norm begins to increase
at some point of the iteration due to the accumulation of the inverted noise,
as shown in Figures 23, 23, and 23. At half precision, the norm grows after
the $11^{th}$ iteration even without noise in Figure 23. The reason this time
is not the attached noise but the truncation errors.

Figure 21: Error norm for size 64 problem with zero noise.

Figure 22: Error norm for size 64 problem with 1% noise.

Figure 23: Error norm for size 64 problem with 10% noise.

##  4 Chebyshev Semi-Iterative Method

###  4.1 Method Overview

The Chebyshev Semi-Iterative (CS) Method requires no inner product
computation, so it is safer in terms of overflow issues when at low precision.
The trade-off here is that the CS method needs prior knowledge of the range of
singular values of the matrix $A$ [14].
As with CG, the CS method is generally derived for symmetric positive definite
linear systems, but there are variations for least squares problems, which are
applied to the normal equations $A^{T}A{\bm{x}}=A^{-1}{\bm{b}}$. The following
algorithm222We are using Algorithm 3 from [11], but we remark that there is a
typographical error for the $\alpha$ parameter for the case $k=1$. We show the
correct formula in our paper. describes the implementation of the CS method
for least squares problems that we used in our work [11] [5]:

Algorithm 2 Chebyshev semi-iterative method

Given $A\in\mathbb{R}^{m\times n}$, ${\bm{b}}\in\mathbb{R}^{m}$, and a
tolerance $\epsilon>0$, choose $0<\sigma_{L}<\sigma_{U}$ such that all nonzero
singular values of A in [$\sigma_{L}$, $\sigma_{U}$], and let
$d=\frac{{\sigma_{U}}^{2}+{\sigma_{L}}^{2}}{2}$ and
$c=\frac{{\sigma_{U}}^{2}-{\sigma_{L}}^{2}}{2}$

Let ${\bm{x}}=0$, ${\bm{v}}=0$, and ${\bf r}={\bm{b}}$.

for $k=0,1,\ldots,\lceil{\log\epsilon-\log
2}/{\log\frac{\sigma_{U}-\sigma_{L}}{\sigma_{U}+\sigma_{L}}}\rceil$ do

$\beta\leftarrow\begin{cases}0,&\text{if }k=0\\\
\frac{1}{2}(c/d)^{2},&\text{if }k=1\\\ (\alpha c/2)^{2},&\text{otherwise,
}\end{cases}\alpha\leftarrow\begin{cases}1/d,&\text{if }k=0\\\
1/(d-c^{2}/(2d)),&\text{if }k=1\\\ 1/(d-\alpha c^{2}/4),&\text{otherwise
}\end{cases}$

${\bm{v}}\leftarrow\beta{\bm{v}}+A^{T}{\bf r}$.

${\bm{x}}\leftarrow{\bm{x}}+\alpha{\bm{v}}$.

${\bf r}\leftarrow{\bf r}-\alpha A{\bm{v}}$.

end for

###  4.2 Experiment

Again, we used image deblurring and tomography reconstruction as two test
problems for the CS method in low precision. At first we ran the image
deblurring test problems of size 32 with default blur and $1\%$ noise in
double precision. We got estimations for the bound of singular values of
matrix $A$ with the built-in MATLAB function SVDS. We set the maximum number
of iterations to be $500$, and plotted the graph at the first, $5^{th}$ (where
the error norm is smallest), and last iteration, as shown in Figures 26, 26,
and 26.

Figure 24: Reconstruction at first iteration.

Figure 25: Reconstruction at the best iteration.

Figure 26: Reconstruction at the last iteration.

The reconstruction is doing poorly even at the best iteration (Figure 26). At
the last iteration, it is even worse due to over-fitting. After a closer look,
we noticed that the number of iterations recommended by the algorithm is at
order of magnitude of $10^{12}$. If we increase the size of the problem to
$64\times 64$, the number overflows to -Inf, which results from the tiny
estimation of the lower bound of matrix $A$’s singular values. The estimation
is so close to zero that we suspect this is actually a result of round-off
errors on a zero entry. Therefore, we implemented Tikhonov regularization for
two reasons: a) to avoid over-fitting; and b) to obtain a more valid lower
bound of singular values by increasing them to a larger value. Specifically,
the singular values for the matrix $A$ after Tikhonov regularization with
regularization parameter $\lambda$ would be
$\sqrt{\sigma_{i}^{2}+\lambda^{2}}$. The technique is discussed in more
details in the next section 4.2.1.

####  4.2.1 Tikhonov Regularization

Tikhonov Regularization includes a regularization term to the original least
squares problem:

$\min_{x}\\{||A{\bm{x}}-{\bm{b}}||_{2}^{2}+\lambda^{2}||{\bm{x}}||_{2}^{2}\\},$

where the regularization parameter $\lambda$ balances the residual term
$||A{\bm{x}}-{\bm{b}}||_{2}^{2}$ and the regularization term
$||{\bm{x}}||_{2}^{2}$. We can rewrite the Tikhonov problem as a least squares
problem

$\min_{{\bm{x}}}\begin{Vmatrix}\begin{pmatrix}A\\\ \lambda
I\end{pmatrix}{\bm{x}}-\begin{pmatrix}{\bm{b}}\\\
0\end{pmatrix}\end{Vmatrix}_{2},$

and the solution to this least squares problem is

${\bm{x}}_{\lambda}=(A^{T}A+\lambda^{2}I)^{-1}A^{T}{\bm{b}}.$

To see why Tikhonov regularization is effective, observe that if we substitute
the singular value decomposition, $\displaystyle A=U\Sigma V^{T}$, into the
expression for ${\bm{x}}_{\lambda}$, and expand the matrix multiplication
column-wise, we would get

${\bm{x}}_{\lambda}=\sum_{i=1}^{n}\phi_{i}^{[\lambda]}\frac{{\bm{u}}_{i}^{T}{\bm{b}}}{\sigma_{i}}{\bm{v}}_{i},$

where the filter factors $\phi_{i}^{[\lambda]}$ are

$\phi_{i}^{[\lambda]}=\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda^{2}}.$

Notice that the filter factor $\phi_{i}^{[\lambda]}$ is approximately equal to
$0$ for smaller singular values and approximately equal to $1$ for larger
singular values. It therefore acts like a filter by decreasing the effects of
magnifying noise in ${\bm{b}}$ when divided by tiny singular values.
With regularization, we are running the CS method with the matrix
$\begin{pmatrix}A\\\ \lambda I\end{pmatrix}$ instead of $A$.
Given $A=U\Sigma V^{T}$, we have

$A^{T}A=V\Sigma^{T}\Sigma V^{T}$

and

$\begin{pmatrix}A^{T}&\lambda I\\\ \end{pmatrix}\begin{pmatrix}A\\\ \lambda
I\end{pmatrix}=A^{T}A+\lambda^{2}I=V\Sigma^{T}\Sigma
V^{T}+V\lambda^{2}V^{T}=V(\Sigma^{T}\Sigma+\lambda^{2})V^{T}.$

Therefore, we know the singular values of $\begin{pmatrix}A\\\ \lambda
I\end{pmatrix}$ are $\sqrt{\sigma_{i}^{2}+\lambda^{2}}$, where $\sigma_{i}$
are singular values of $A$. When $\sigma_{i}$ is small,
$\sqrt{\sigma_{i}^{2}+\lambda^{2}}\approx\lambda$, and we can directly use
$\lambda$ as the lower bound of singular values for the matrix after
regularization. With Tikhonov regularization, the singular values are
increased to a more reasonable lower bound.

####  4.2.2 Image Deblurring

Figure 27: Image reconstruction in double precision of size 64 problem with
mild blurring and no noise.

Figure 28: Image reconstruction in single precision of size 64 problem with
mild blurring and no noise.

Figure 29: Image reconstruction in half precision of size 64 problem with mild
blurring and no noise.

With regularization, we successfully ran $478$ iterations in all three
precisions, without the occurrence of NaNs. We plotted the results in Figures
29, 29, and 29. For half precision, the result is not clear as expected, and
the image seems to be dominated by round-off errors. We plotted the error norm
in Figure 35.

The error norms do not converge to a point like those generated by the CGLS
algorithm. Instead, they oscillate at the beginning, and then converge for
double and single precision. However, for half precision, the error norm
increases rapidly, indicating that the round-off errors accumulate and take
over. We suspect this is because the regularization parameter is too small, so
we may need to develop better ways to determine a suitable regularization
parameter for half precision.

We then added noise to ${\bm{b}}$ and displayed the resulting computed
reconstructions in Figures 32, 32, and 32.

Figure 30: Image reconstruction in double precision of size 64 problem with
mild blurring and 10% noise.

Figure 31: Image reconstruction in single precision of size 64 problem with
mild blurring and 10% noise.

Figure 32: Image reconstruction in half precision of size 64 problem with mild
blurring and 10% noise.

Here we only showed results with $10\%$ noise; results for other noise levels
are consistent with those shown in previous sections. Surprisingly, the result
in half precision looks similar to the one in double or single precision, and
it is much better than the image without noise. Moreover, the clear result is
acquired by only $16$ iterations, illustrating the efficiency of the CS
method.

Figure 33: Error norm of a size 64 problem with mild blurring of different
precisions with 0 noise.

Figure 34: Error norm of a size 64 problem with mild blurring of different
precisions with 1% noise.

Figure 35: Error norm of a size 64 problem with mild blurring of different
precisions with 10% noise.

By observing the error norm plots in Figure 35, we also see that there seems
to be very little difference in the convergence behavior of CS for the
different precision levels when the noise level is high. The more noise we
have, the closer are the results between double and half precision.

####  4.2.3 Tomography Reconstruction

Although CS avoids the computation of inner products, overflow still occurs at
half precision during the calculation of matrix-vector multiplications for the
tomography reconstruction problem. NaNs start to occur at the second
iteration. Figure 48 shows the resulting image at the end of the first
iteration at half precision for a problem with $10\%$ noise, which is already
reasonably good with visible shapes and boundaries.

Figure 36: Half precision, size 64, 10% noise, without rescale.

Due to overflow, we rescaled the problem as for CGLS so that CS could run more
iterations for half precision. The algorithm successfully runs to the end
without occurrences of overflow. In the figures below, we show the results at
the last iteration for problems with different precision and noise levels.

Figure 37: Double precision, size 64, zero noise.

Figure 38: Double precision, size 64, 1% noise.

Figure 39: Double precision, size 64, 10% noise.

Figure 40: Single precision, size 64, zero noise.

Figure 41: Single precision, size 64, 1% noise.

Figure 42: Single precision, size 64, 10% noise.

Figure 43: Half precision, size 64, zero noise (after rescaling).

Figure 44: Half precision, size 64, 1% noise (after rescaling).

Figure 45: Half precision, size 64, 10% noise (after rescaling).

After rescaling, when there is noise on the right hand side ${\bm{b}}$ of the
problem, the resulting image at low precision is as valid as images produced
at high precision. With $10\%$ noise (Figures 45, 45, and 45), the CS
algorithm ran $32$ iterations, and produced good reconstructions with clear
boundaries and background for all three precision levels.

However, when there is zero noise, the result at half precision is a poor
reconstruction. The background and the object both look blurry and blend
together, and the object contents are hardly visible. We believe the
dissatisfying image is a combined result of accumulation of round-off errors
and under-regularization.

After a closer look at the result of each iteration, we noticed that the
reconstructions in the first few iterations look smooth and improved as the
iteration moves on, but at some point they start to become noisy and blurry.
Unlike CGLS where the output image at each iteration refines steadily, the
results from CS have more oscillations from iteration to iteration, and it has
a general trend of becoming noisy as the iteration goes on. Below we show the
resulting images at the first iteration (Figure 48), the $5^{th}$ iteration
(Figure 48) where the error norm is the smallest, and the $250^{th}$ iteration
(Figure 48). At the $250^{th}$ iteration the image is already very noisy as
round-off errors accumulated along the way.

Figure 46: Reconstruction at first iteration, half precision, zero noise
(after rescaling).

Figure 47: Reconstruction at best iteration, half precision, zero noise (after
rescaling).

Figure 48: Reconstruction at $250^{th}$ iteration, half precision, zero noise
(after rescaling).

Figures 51, 51 and 51 present the error norms at different precision levels
across different noise levels. As expected, the error norms overlap for test
problems with noise, implying that the quality of the reconstructions is
similar. However, for the noise-free problems, the error norms follow a
decreasing trend for double and single precision and a slightly increasing
trend for half precision. Besides, the norms oscillate more as noise level
decreases, which corresponds with our observation from the resulting images.

Figure 49: Error norm of a size 64 problem at different precisions with zero
noise.

Figure 50: Error norm of a size 64 problem at different precisions with 1%
noise.

Figure 51: Error norm of a size 64 problem at different precisions with 10%
noise.

##  5 Discussion

We incorporated CGLS with regularization to compare its performance with that
of CS more fairly. The effect of regularization is apparent, as the resulting
image is less noisy, especially for low precision. As the noise level
increases, the difference between high precision and low precision becomes
less significant. Yet still, rescaling is necessary for half precision to
avoid overflow for the tomography reconstruction problem. Below we showed the
result of CGLS with regularization for the two test problems.

Figure 52: CGLS with regularization, double precision, 10% noise.

Figure 53: CGLS with regularization, single precision, 10% noise.

Figure 54: CGLS with regularization, half precision, 10% noise.

Figure 55: CGLS with regularization, double precision, zero noise.

Figure 56: CGLS with regularization, single precision, zero noise.

Figure 57: CGLS with regularization, half precision, zero noise (after
rescaling).

Figure 58: CGLS with regularization, double precision,10% noise.

Figure 59: CGLS with regularization, single precision,10% noise.

Figure 60: CGLS with regularization, half precision,10% noise (after
rescaling).

Overall, CGLS is a more stable algorithm than CS and has less oscillations in
the solution. The result given by CGLS is steadily improving in the iterations
even with the presence of noise. Switching to lower precision does not have a
large impact on the performance of CGLS, as long as no overflow or underflow
occurs. And if they do occur, we can rescale the problem to delay the
occurrence of overflow and NaNs so that a relatively good result can still be
obtained. However, we do expect there to be still significant problems with
overflow and NaNs for larger problems.

The main problem with CGLS at low precision is to find the suitable rescaling
parameter, which depends on the matrix $A$ and vector ${\bm{b}}$ of each
problem. We did not have issues of overflow/underflow at all for the image
deblurring problem, while NaNs started to appear at the first iteration for
the tomography reconstruction problem. We tried several rescaling parameter
before we found a suitable one, and sometimes it is hard to determine whether
a parameter is suitable or not if we do not have an idea of what the real
image looks like. The purpose of rescaling is to avoid NaNs in the solution,
but a solution without NaNs does not imply it has been rescaled properly. For
example, Figure 61 is the result of a size-$32$ tomography reconstruction test
problem given by CGLS after $100$ iterations. No NaN occurred during the
iteration process, but the resulting image is far from the true solution. This
is because the overflow leads to an underflow of ${\bm{x}}$ to zero
immediately in the first iteration. Though it did not result in NaNs, the
algorithm did not capture any information about ${\bm{x}}$ in that iteration
as well. Therefore, the image is not rescaled properly despite the absence of
NaNs.

Figure 61: Half precision, size 32, zero noise.

CS, on the other hand, does not necessarily require rescaling at low precision
as it avoids the calculation of inner products, which is the main source of
overflow. The major advantage of CS is that it is able to run more iterations
before the first NaNs occur, so that the resulting images have more chances to
get refined and are less likely to be cut off before the optimal point. There
are also cases when CGLS overflows at the first iteration and produces
anything but NaNs, while CS is able to run several iterations and produce a
relatively meaningful image, as for the tomography reconstruction test
problem. In these circumstances, if only a rough estimate of the original
image is desired, CS can be a more convenient choice as it saves the trouble
of finding a proper rescaling parameter.

However, CS requires the user to have an estimate of the bounds of $A$’s
singular values, which may be hard to estimate. Even with Tikhonov
regularization, the resulting bound depends on the regularization parameter
$\lambda$, which in turn depends on the noise level and the problem itself.
When $\lambda$ is small, the computation becomes more risky, especially at low
precision. Generally speaking, when the noise is large, we would need a larger
$\lambda$ to filter out the noise. Therefore, it is not surprising that CS has
better performance when the noise level is high. However, when the noise is
negligible, the more stable CGLS is a better choice.

##  6 Conclusion

In this project, we explored methods for solving inverse problems at low
precision. We first modified several MATLAB built-in operations using chop for
custom precision levels and applied the blocking method to decrease the error
bounds of the simulation. Then we ran the modified CGLS and CS methods on
image deblurring and tomography reconstruction tasks to compare their
performance across different precision levels. We added Tikhonov
regularization to both methods to balance signal and noise.

After comparing results given by the two algorithms, we concluded that CGLS is
more stable than CS. Its output image steadily gets closer to the real image.
However, CGLS is likely to suffer from overflow at low precision as it
involves the calculation of inner products. One solution to this issue is to
rescale the problem, but finding a suitable rescaling parameter requires trial
and error. Moreover, the rescaling is unlikely to be effective for extremely
large problems when using low precision. In the future, we hope to find a more
direct way that finds a suitable rescaling parameter given the matrix $A$,
vector ${\bm{b}}$, and precision level. Based on intuition gained from
experiments, a factor that “rescales” the values to one digit is often a good
choice. We would also look at algorithms for normalizing $A$ as described in
the survey by Higham and Mary [9].

The performance of the CS method is less stable and depends more on the noise
level of the problem. It has more oscillations than CGLS during the iterating
process. However, CS requires no calculation of inner products, and it is
therefore less likely to overflow. The method avoids the need of rescaling for
some problems; with rescaling, it can refine the result with even more
iterations and often ends up with a better resulting image than CGLS. The
price paid here is that CS needs some prior knowledge of the range of the
matrix’s singular values. When the matrix is large and sparse, for low
precision levels, we need a regularization parameter large enough to obtain a
valid lower bound and avoid the accumulation of round-off errors. For problems
with noise on the right hand side ${\bm{b}}$, most of the time we would
naturally end up with a large enough $\lambda$. Therefore, CS has great
performance for noisy problems. However, when the problem is noise-free, the
estimated $\lambda$ given by the parameter selection methods are often too
small, and CS performs poorly. In the future, we hope to develop better
methods for choosing suitable regularization parameters for CS at low
precision.

When examining the error norms of the solution at each iteration, we noticed
that sometimes the error norm did not match the so-called “eyeball norm.” An
image with clear shape and background may have a higher error norm than a
blurry, noisy image. The latter one, though with a smaller error norm, is
obviously less informative than the former one. Therefore, more research is
needed on ways to take other aspects of the output image into the account of
error measurement so that we could find an image that conveys the most
information.

Moreover, we would like to explore more iterative methods other than CGLS and
CS and implement them in low precision, as well as mixed precision, which is
expected to have both the accuracy of high precision and the gains of speed
from low precision.

## References

  * [1] Å. Björck.  Numerical Methods for Least Squares Problems.  SIAM, 1996.
  * [2] C. L. Epstein.  Introduction to the Mathematics of Medical Imaging, Second Edition.  SIAM, Philadelphia, PA, 2007.
  * [3] S. Gazzola, P. C. Hansen, and J. G. Nagy.  IR tools: a MATLAB package of iterative regularization methods and large-scale test problems.  Numerical Algorithms, 81(3):773–811, 2019.
  * [4] G. Gupta.  Using tensor cores for mixed-precision scientific computing.  NVIDIA Technical Blog, https://developer.nvidia.com/blog/tensor-cores-mixed-precision-scientific-computing/, Oct-2021.
  * [5] M. H. Gutknecht and S. Röllin.  The Chebyshev iteration revisited.  Parallel Computing, 28(2):263–283, 2002.
  * [6] P. C. Hansen.  Discrete Inverse Problems: Insight and Algorithms.  SIAM, 2010.
  * [7] M. R. Hestenes and E. Stiefel.  Methods of conjugate gradients for solving.  Journal of research of the National Bureau of Standards, 49(6):409, 1952.
  * [8] N. J. Higham.  Accuracy and Stability of Numerical Algorithms.  SIAM, 2002.
  * [9] N. J. Higham and T. Mary.  Mixed precision algorithms in numerical linear algebra.  Acta Numerica, 31:347–414, 2022.
  * [10] N. J. Higham and S. Pranesh.  Simulating low precision floating-point arithmetic.  SIAM Journal on Scientific Computing, 41(5):C585–C602, 2019.
  * [11] X. Meng, M. A. Saunders, and M. W. Mahoney.  LSRN: A parallel iterative solver for strongly over-or underdetermined systems.  SIAM Journal on Scientific Computing, 36(2):C95–C118, 2014.
  * [12] P. San Juan, R. Rodríguez-Sánchez, F. D. Igual, P. Alonso-Jordá, and E. S. Quintana-Ortí.  Low precision matrix multiplication for efficient deep learning in nvidia carmel processors.  The Journal of Supercomputing, 77(10):11257–11269, 2021.
  * [13] L. N. Trefethen and D. Bau III.  Numerical Linear Algebra, volume 50.  SIAM, 1997.
  * [14] H. Wang.  A Chebyshev semi-iterative approach for accelerating projective and position-based dynamics.  ACM Transactions on Graphics (TOG), 34(6):1–9, 2015.


there is no path from $s$ to $t$ then terminates.

    2. (b)

Construct a block flow $f$ in $G_{L}$ using Edmonds–Karp algorithm.

    3. (c)

Do 3a and 3b for the residual network $G_{f}$.

######  Theorem 24 (Block flow).

Consider a block flow $f$. Let’s denote the minimum number of edges from $s$
to a vertex $v$ in $G$ and in $G_{f}$ by $dist(s,v)$ and $dist^{\prime}(s,v)$
respectively. Then, $dist^{\prime}(s,t)\geq dist(s,t)+1$.

[1em]0em

###### Proof.

Let’s denote by $len(P)$ the length of a path $P$. First let’s note that if
there is a path $P^{\prime}$ in $G_{f}$ from $s$ to $t$ then either
$len(P^{\prime})>dist(s,t)$ or it contains a new edge $e\notin E(G)$. Thus it
is sufficient to prove that for any path $P^{\prime}$ in $G_{f}$ from $s$ to
$t$ that contains a new edge $e\notin E(G)$ the length $len(P^{\prime})\geq
dist(s,t)+1$.

Let’s proof moreover: if the path $P^{\prime}$ in $G_{f}$ from $s$ to any
vertex $v$ contains new edges $e_{1},e_{2},...,e_{k}\notin E(G)$ then
$len(P^{\prime})\geq dist(s,v)+2k$. Let’s prove it by induction by $k$:

    1. (a)

If the path $P^{\prime}$ in $G_{f}$ from $s$ to a vertex $v$ contains a new
edge $e_{1}=(w_{1},u_{1})$ thus the inverse edge $(u_{1},w_{1})$ belongs to
$G$ and $f(u_{1},w_{1})=c(u_{1},w_{1})$ for the block flow $f$. By
decomposition algorithm and previous statement there exists a path in the
level graph $G_{L}$ that contains this edge $(u_{1},w_{1})$. Thus
$dist(s,w_{1})=dist(s,u_{1})+1$.

Let’s denote the part of the path $P^{\prime}$ from the vertex $a$ to $b$ by
$P^{\prime}_{ab}$. Since $P^{\prime}_{sw_{1}}$ and $P^{\prime}_{u_{1}v}$ are
paths in $G$ thus $len(P^{\prime}_{sw_{1}})\geq dist(s,w_{1})$ and
$len(P^{\prime}_{u_{1}v})\geq dist(u_{1},v)$. Thus
$len(P^{\prime}_{su_{1}})\geq dist(u_{1})+2\implies len(P^{\prime})\geq
dist(s,v)+2$.

    2. (b)

Let this statement holds for any $j\leq k-1$ then let’s prove it for $j=k$.
Since $e_{k}=(w_{k},u_{k})$ is a new edge thus there exists a path in the
level graph $G_{L}$ that contains the inverse edge $(u_{k},w_{k})$ and thus
$dist(s,w_{k})=dist(s,u_{k})+1$. Since $P^{\prime}_{u_{k}v}$ is a path in $G$
thus $len(P^{\prime}_{u_{k}v})\geq dist(u_{k},v)$.

By the induction statement the length $len(P^{\prime}_{sw_{k}})\geq
dist(s,w_{k})+2(k-1)$. Thus $len(P^{\prime}_{su_{k}})\geq
dist(s,u_{k})+2k\implies len(P^{\prime})\geq dist(s,v)+2k$.

Hence, if a path $P^{\prime}$ in $G_{f}$ from $s$ to $t$ contains a new edge
$e\notin E(G)$ then the length $len(P^{\prime})\geq dist(s,t)+1$. ∎

NB! By the block flow theorem 24 the minimum path length in Edmonds–Karp
algorithm is increased at least by 1 at each step and thus the number of steps
is $O(\|V\|)$. By the decomposition theorem 22 the number of steps in
Edmonds–Karp algorithm for a level graph $G_{L}$ is $O(\|V\|\|E\|)$ (like in
decomposition algorithm). Thus the complexity of Dinic algorithm is
$\mathbf{O(\|V\|^{2}\|E\|)}$.

NB! One can improve the complexity of Dinic algorithm using the dynamic tree
structure and in that case it will be $\mathbf{O(\|V\|\|E\|\log(\|V\|))}$.

###  14.3 Proof of Menger’s theorem

######  Lemma 42 (Integer flow).

If all capacities of a flow network are integer then there exists a max-flow
that is also integer on every edge.

[1em]0em

###### Proof.

By using any max-flow construction algorithm one can obtain the integer max-
flow. ∎

######  Lemma 43.

If all capacities of a flow network are either 0 or 1 and the max-flow value
equal to $K$ then the paths from the source $s$ to the terminal $t$ are
$K$-edge independent.

[1em]0em

###### Proof.

Let’s construct the integer max-flow from the previous lemma. By using
decomposition algorithm any paths $P_{i}$ will be edge independent and since
the flow among each path equals to 1 then the number of paths $P_{i}$ is $K$.
∎

######  Theorem 25 (Menger).

A graph is $\kappa$-vertex (edge) connected _iff._ for any two vertices there
exist $\kappa$-vertex (edge) independent paths.

[1em]0em

###### Proof.

  1. 1.

Let’s prove first the edges condition.

    1. (a)

_Sufficiency._ If there exist $\kappa$-edge independent paths thus there is no
$\kappa-1$-edge cuts.

    2. (b)

_Necessarity._ Let the graph $G$ be connected after a deletion of any
$\kappa-1$ edges. Let’s denote two vertices in the condition by $s$ and $t$.
Let’s define weight of every edge in $G$ by 1. Consider $G$ as the flow
network and let $f$ be the max-flow. By the Ford-Fulkerson theorem 23:
$F=c(A_{min},B_{min})$ and equal to the number of edges from vertices of
$A_{min}$ to vertices of $B_{min}$. Since the graph $G$ be connected after a
deletion of any $\kappa-1$ edges and $F$ is integer then $F\geq\kappa$. Thus
by the previous lemma the paths from $s$ to $t$ are $\kappa$-edge independent.

  2. 2.

Let’s proof the vertex condition.

    1. (a)

_Sufficiency._ The same as for edges.

    2. (b)

For _necessarity_ condition let’s give additional definition:

######  Definition 14.11.

Let $A$ and $B$ be two subsets of vertices of the graph $G$. A set $S\subset
V(G)$ is called $AB$-separator if there is no path in $G$ starting from any
vertex of $A$ and ending in any vertex of $B$ after the deletion of the set
$S$.

Let’s proof that for any $A\subset V(G)$ and $B\subset V(G)$ if any
$AB$-separator consist of at least $\kappa$ vertices then there exist
$\kappa$-vertex independent paths from vertices of $A$ to vertices $B$ (here
paths can consist of just one vertex). Let’s proof it by induction of number
of edges in $G$:

      1. i.

If $G$ has no edges then the vertices in $A\bigcap B$ are also paths.

      2. ii.

Let we know the induction statement for the graph
$G^{\prime}=G\setminus(v_{1},v_{2})$. Let’s proof it for $G$. Let’s assume
that $AB$-separator in $G$ consists of $\kappa$ vertices. If $AB$-separator in
$G^{\prime}$ consists also of $\kappa$ vertices then by the induction
statement there exist $\kappa$-vertex independent paths. Thus there exists
$AB$-separator $S$ in $G^{\prime}$ that consists of $\kappa-1$ vertices. Let
$v_{1}$ (without loss of generality) is reachable from some vertex of $A$ in
$G^{\prime}\setminus S$ (otherwise $S$ is a $AB$-separator in $G$).

Let’s denote this set $S\bigcup\\{v_{1}\\}$ by $S_{1}$. The set $S_{1}$ is
$AB$-separator in $G$ (otherwise $S$ is not a $AB$-separator in $G^{\prime}$).
The $AS_{1}$-separator is also $AB$-separator in $G$ and thus consists of at
least $\kappa$ vertices thus by the induction statement there exist
$\kappa$-vertex independent paths from vertices of $A$ to vertices $S_{1}$. By
providing the same proof for $S_{2}=S\bigcup\\{v_{2}\\}$ one can obtain
$\kappa$-vertex independent paths in $G$ from vertices of $S_{2}$ to $B$.
Since $\|S_{1}\|=\|S_{2}\|=\kappa$ there exist $\kappa$-vertex independent
paths from $A$ to $B$ in $G$.

∎

###  14.4 Minimum-cost flow

######  Definition 14.12.

Let a cost function $a:V(G)\times V(G)\rightarrow\mathbb{R}$ and a flow $f$ is
defined on a flow network $G$. The total cost of the flow is
$p(G)=\sum\limits_{(u,v)\in V(G)\times V(G),\ f(u,v)>0}a(u,v)f(u,v)$.

######  Definition 14.13.

A flow is called min-cost-flow if $f$ has the minimal total cost $p$ among the
all another flows with the flow value $F$.

Let $G$ be a flow network with a flow $f$ and cost function $a$.

######  Lemma 44 (Necessary and sufficient condition of min-cost-flow).

A flow $f$ is min-cost-flow iff. there are no negative cycles corresponding to
the cost function $a$ in the residual network $G_{f}$.

[1em]0em

###### Proof.

  1. 1.

_Necessarity._ Assume the contrary: there exists negative cycle $C$ in the
residual network $G_{f}$. Let $c_{f,C}=\min_{e\in C}c_{f}(e)$. Let’s define a
flow

$f_{C}(u,v)=\begin{cases}c_{f,C}&(u,v)\in C\\\ -c_{f,C}&(v,u)\in C\\\
0&otherwise\end{cases}.$

Thus the total cost of $f+f_{C}$ will be less than the total cost of $f$.

  2. 2.

_Sufficiency._ Let $f$ be the flow with no negative cycles corresponding to
the cost function $a$ in the residual network $G_{f}$ and $f^{\prime}$ min-
cost-flow among the flows with flow value $F$. Thus $p(f^{\prime})\leq p(f)$.
By the subtraction of flows lemma 38:
$f^{\prime}=f+\sum\limits_{i}f_{C_{i}}\implies p(f^{\prime})\geq p(f)\implies
p(f^{\prime})=p(f)$.

∎

######  Theorem 26.

Let $f$ be the min-cost-flow, $P$ be the path with minimal sum of costs among
the paths from the source $s$ to the terminal $t$ in the residual network
$G_{f}$ and $c_{f,P}$ be the minimal residual capacity among the path $P$. Let

$\forall\delta:0\leq\delta\leq c_{f,P}\quad
f_{\delta}(u,v)=\begin{cases}\delta&(u,v)\in P\\\ -\delta&(v,u)\in P\\\
0&otherwise\end{cases}.$

Then $f+f_{\delta}$ is a min-cost-flow among the flows with value $F+\delta$.

[1em]0em

###### Proof.

Let $g$ be a min-cost-flow among the flows with value $F+\delta$. Then by the
two flow lemma 37: $g=f+f^{\prime}$ there $f^{\prime}$ is a flow in the
residual network $G_{f}$ with a flow value $\delta$. By the decomposition
theorem 22
$f^{\prime}=\sum\limits_{i}f^{\prime}_{P_{i}}+\sum\limits_{i}f^{\prime}_{C_{i}}$.

By previous lemma any cycle $C_{i}$ can not be negative cycle corresponding to
the cost function $a$. If $C_{i}$ is positive cycle then one can decreased the
total cost of $g$ and it holds contradiction. Thus $\forall i,e\in
C_{i}:a(e)=0$.

$p(f^{\prime})=\sum\limits_{i}p(f^{\prime}_{P_{i}})=\sum\limits_{i}\sum\limits_{(u,v)\in
P_{i}}a(u,v)f^{\prime}_{P_{i}}(u,v)=\sum\limits_{i}f^{\prime}(P_{i})\sum\limits_{(u,v)\in
P_{i}}a(u,v)\geq$ $\geq\sum\limits_{i}f^{\prime}(P_{i})\sum\limits_{(u,v)\in
P}a(u,v)=\delta\sum\limits_{(u,v)\in P}a(u,v)=p(f_{\delta}).$

∎

Before the considering min-cost-max-flow construction algorithm let’s define
the flow among a path $P$ in the residual network by

$f^{\prime}=\begin{cases}c_{f,P}&(u,v)\in P\\\ -c_{f,P}&(v,u)\in P\\\
0&otherwise\end{cases},$

where $c_{f,P}$ is the minimal residual capacity among the path $P$.

  1. 1.

Min-cost-max-flow construction using Ford-Fulkerson algorithm.
_Description:_

    1. (a)

Initialise $f=0$.

    2. (b)

Find a negative cycle $C$ corresponding to the cost function $a$ in the
residual network $G_{f}$ using Belman-Ford algorithm.

    3. (c)

Increase $f$ by the flow $f^{\prime}$ among the cycle $C$ and go to step 1b.

    4. (d)

Find the path $P$ with minimal sum of costs from the source $s$ to the
terminal $t$ in the residual network $G_{f}$.

    5. (e)

Increase $f$ by the flow $f^{\prime}$ among the path $P$ and go to step 1d.

NB! Since this algorithm is modification of Ford-Fulkerson algorithm the
complexity is also depends of the distribution of capacities and price
function.

  2. 2.

Min-cost-max-flow construction using Dinic algorithm.
_Description:_

    1. (a)

Find max-flow $f$ using Dinic algorithm.

    2. (b)

Find a negative cycle $C$ corresponding to the cost function $a$ in the
residual network $G_{f}$ using Belman-Ford algorithm.

    3. (c)

Increase $f$ by the flow $f^{\prime}$ among the cycle $C$ and go to step 2b.

NB! The complexity of this algorithm is
$\mathbf{O(\|V\|\|E\|^{2}+\|V\|^{2}\|E\|)=O(\|V\|\|E\|(\|E\|+\|V\|))}.$

##  15 Local and global characteristics of graph

Consider some local characteristics of an undirected graph $G$.

######  Definition 15.1.

Centrality is a function defined for each vertex of a graph that contains some
information of a graph structure.

First let’s give some general centralities which contain local information of
the graph $G$:

  1. 1.

Let’s denote by $N(v)$ the set of vertices which adjacent to the vertex $v$.
The simplest example is the degree centrality $\deg{v}=\|N(v)\|$.

  2. 2.

Let’s denote by $G(N(v))$ the maximal sub-graph on vertices $V(N(v))$. Then
$MC(v)$ be the largest connected component in $G(N(v))$. Maximum neighborhood
component $MNC(v)$ is the number of vertices in $MC(v)$.

  3. 3.

Density of maximum neighborhood component
$DMNC(v)=\frac{\|E(MC(v))\|}{\|V(MC(v))\|^{\epsilon}}$, for some
$\epsilon\in[1,2]$.

  4. 4.

Local cluster coefficient $C(v)=\frac{\text{number of edges in
}N(v)}{\text{maximum possible number of edges in
}N(v)}=\frac{2\|E(N(v)))\|}{\deg(v)(\deg(v)-1)}$.

Next the list of general global characteristics of a graph is represented:

  1. 1.

The simplest example is the diameter $diam(G)=\max_{s,t}dist(s,t)$.

  2. 2.

Density $D(G)=\frac{\text{number of edges in }G}{\text{maximum possible number
of edges in }G}=\frac{2\|E(G))\|}{\|V(G)\|(\|V(G)\|-1)}$.

  3. 3.

Global efficiency
$E_{glob}(G)=\frac{1}{\|V(G)\|(\|V(G)\|-1)}\sum\limits_{s\neq
t}\frac{1}{dist(s,t)}$.

  4. 4.

Average shortest path length
$L(G)=\frac{1}{\|V(G)\|(\|V(G)\|-1)}\sum\limits_{s\neq t}dist(s,t)$.

  5. 5.

Average clustering coefficient $C(G)=\frac{1}{\|V(G)\|}\sum\limits_{v\in
V(G)}C(v)=\frac{1}{\|V(G)\|}\sum\limits_{v\in
V(G)}\frac{2\|E(N(v)))\|}{\deg(v)(\deg(v)-1)}$.

  6. 6.

Small world coefficient
$SW(G)=\frac{C(G)}{C(G_{rand})}/\frac{L(G)}{L(G_{rand})}$, where $G_{rand}$ is
a random graph $(\|V(G)\|,\|E(G)\|)$.

NB! It turns out that most real life networks are satisfied the small world
property $SW(G)\gg 1$.

Here the list of general centralities contained a global information of a
graph $G$ is represented:

  1. 1.

Betweenness centrality $BC(v)=\sum\limits_{s\neq t\neq
v}\frac{\sigma_{st}(v)}{\sigma_{st}}$, where $\sigma_{st}$ is the total number
of shortest paths from $s$ to $t$ and $\sigma_{st}(v)$ is the total number of
shortest paths which contains vertex $v$.

  2. 2.

Let $T_{s}$ be a shortest path tree with the root $s$. Bottleneck
$BN(v)=\sum\limits_{s\in V(G)}p_{s}(v)$, where

$p_{s}(v)=\begin{cases}1&\text{if the set of descendants of $v$ in $T_{s}$
contains more than $\|V(T_{s})\|/4$ vertices}\\\
0&\text{otherwise}\end{cases}.$

  3. 3.

Closeness centrality $Clo(v)=\frac{1}{\sum\limits_{t\in V(G)}dist(v,t)}$.

  4. 4.

Eccentricity centrality $EC(v)=\frac{1}{\max_{t\in V(G)}dist(v,t)}$.

  5. 5.

Local efficiency $E_{loc}(G)=\frac{1}{\|V(G)\|}\sum\limits_{v\in
V(G)}E_{glob}(N(v)).$

  6. 6.

Radiality $Rad(v)=\frac{\sum\limits_{t\in
V(G)}(diam(G)+1-dist(v,t))}{\|V(G)\|-1}$.

  7. 7.

Stress $Str(v)=\sum\limits_{s\neq t\neq v}\sigma_{st}(v)$, where $\sigma_{st}$
is the total number of shortest paths from $s$ to $t$ which contains vertex
$v$.

First let’s note that if $A(G)$ is the adjacency matrix of $G$ and $E(G)$ is
the efficiency matrix, where

$e_{ij}=\begin{cases}\frac{1}{dist(v_{i},v_{j})}&i\neq j\\\
0&\text{otherwise}\end{cases},$

then

$C(v)=\frac{\sum\limits_{(v_{i},v_{j})\in
E(N(v))}a_{i,j}}{\|V(N(v))\|(\|V(N(v))\|-1)},$
$E_{glob}(N(v))=\frac{\sum\limits_{(v_{i},v_{j})\in
E(N(v))}e_{i,j}}{\|V(N(v))\|(\|V(N(v))\|-1)}.$

Thus

$A(G)-E(G)=\begin{cases}\frac{1}{dist(v_{i},v_{j})}&(v_{i},v_{j})\notin
E(G)\\\ 0&\text{otherwise}\end{cases}.$

Let’s denote by $k$ the number of vertices in $N(v)$. Now let’s prove some
relationship between this measures corresponding to [7]:

######  Theorem 27.

$E_{loc}(G)\leq\frac{1}{2}(1+C(G))$.

[1em]0em

###### Proof.

Let’s prove it for a local characteristics:

$C(v)=\frac{\sum\limits_{(v_{i},v_{j})\in E(N(v))}1}{k(k-1)},$

$\frac{1}{2}(1+C(v))=\frac{1}{2}\frac{\sum\limits_{(v_{i},v_{j})\in
E(N(v))}1+\sum\limits_{(v_{i},v_{j})\in
E(N(v))}1+\sum\limits_{(v_{i},v_{j})\notin
E(N(v))}1}{k(k-1)}=\frac{\sum\limits_{(v_{i},v_{j})\in
E(N(v))}1+\sum\limits_{(v_{i},v_{j})\notin E(N(v))}\frac{1}{2}}{k(k-1)},$

$E_{glob}(N(v))=\frac{1}{k(k-1)}\sum\limits_{s\neq t\in
V(N(v))}\frac{1}{dist(s,t)}\leq\frac{1}{k(k-1)}(\sum\limits_{(v_{i},v_{j})\in
E(N(v))}1+\sum\limits_{(v_{i},v_{j})\notin
E(N(v))}\frac{1}{2})=\frac{1}{2}(1+C(v)).$

∎

######  Theorem 28.

$L(G)\geq 2-D(G).$

[1em]0em

###### Proof.

Let $\alpha$ be the number of pairs of vertices $v,w$ with $dist(v,w)=2$ and
$\beta$ be the number of pairs of vertices $v,w$ with $dist(v,w)\geq 3$ and
let’s denote the number of vertices of $G$ by $n$ and the number of edges by
$m$. Thus $2m+\alpha+\beta=n(n-1)$.

$D(G)=\frac{2m}{n(n-1)},$ $L(G)=\frac{2m+2\alpha+t\beta}{n(n-1)},\text{ for
some }t\geq 3,$
$\frac{2m+2\alpha+t\beta}{n(n-1)}\geq\frac{2m+2(\alpha+\beta)}{n(n-1)}=\frac{4m+2(\alpha+\beta)}{n(n-1)}-\frac{2m}{n(n-1)}=2-D(G).$

∎

######  Theorem 29.

$3-L(G)\leq 2E_{glob}(G)\leq 1+D(G).$

[1em]0em

###### Proof.

Let’s denotations be the same as in previous theorem.

  1. 1.

Case: $2E_{glob}(G)\leq 1+D(G).$

$2E_{glob}(G)=\frac{1}{n(n-1)}(4m+\alpha+2\epsilon\beta),\text{ for some
}\epsilon<\frac{1}{2},$
$\frac{1}{n(n-1)}(4m+\alpha+2\epsilon\beta)\leq\frac{1}{n(n-1)}(4m+\alpha+\beta)=\frac{2m}{n(n-1)}+\frac{2m+\alpha+\beta}{n(n-1)}=1+D(G).$

  2. 2.

Case: $2E_{glob}(G)\geq 3-L(G).$

$\frac{1}{n(n-1)}(4m+\alpha+2\epsilon\beta)\vee
3-\frac{2m+2\alpha+t\beta}{n(n-1)},$ $4m+\alpha+2\epsilon\beta\vee
3(2m+\alpha+\beta)-2m-2\alpha-t\beta,$ $\beta(2\epsilon+t-3)\geq 0.$

∎

NB! The equality in the previous theorems 27-29 is obtained in the case then
$diam(N(v))\leq 2$ for any vertex $v$.

##  16 Laplacian spectrum

###  16.1 The normalized Laplacian and its basic spectral property

Let us introduce the normalized Laplacian on a undirected unweighted graph $G$
with no self-loops.

######  Definition 16.1.

The normalized Laplacian matrix $\Delta$ is defined as
$D^{-\frac{1}{2}}LD^{-\frac{1}{2}}$, where $L$ is the Laplacian matrix of $G$,
and $D=\mathrm{diag}(\deg v_{1},\cdots,\deg v_{n})$ is the diagonal matrix
consisting of degrees of $G$. Then the $(i,j)$-entry of $\Delta$ is

$\begin{cases}1,&\text{ if }i=j\\\ -\frac{1}{\sqrt{d_{i}d_{j}}},&\text{ if
}i\neq j,\\{i,j\\}\in E(G)\\\ 0,&\text{ otherwise},\\\ \end{cases}$

in which we simply write $d_{i}=\deg v_{i}$ for convenience.

We now list the important properties of $\Delta$:

######  Lemma 45 (Properties of normalized Laplacian matrix).

  1. 1.

By viewing $\Delta$ as a linear operator
$\Delta:\mathbb{R}^{n}\to\mathbb{R}^{n}$, $\Delta$ is self-adjoint with
respect to the scalar product $(\cdot,\cdot)$, i.e.,

$(x,\Delta y)=(\Delta x,y)$

for all $x,y\in\mathbb{R}^{n}$. Here the scalar product is defined for pairs
of vectors, formally, for any $x,y\in\mathbb{R}^{n}$,
$(x,y):=\sum_{i=1}^{n}x_{i}y_{i}$.

  2. 2.

$\Delta$ is non-negative:

$(\Delta x,x)\geq 0$

for all $x$.

  3. 3.

$\Delta x=0$ precisely when $x$ is a vector colinear to
$(\sqrt{d_{1}},\cdots,\sqrt{d_{n}})$.

  4. 4.

The trace of $\Delta$ is $n$.

[1em]0em

###### Proof.

  1. 1.

Since $L$ is a symmetric matrix, the normalized Laplacian matrix
$\Delta:=D^{-\frac{1}{2}}LD^{-\frac{1}{2}}$ is also symmetric. Thus, by
regarding a vector $x$ as a column vector,

$(x,\Delta y)=x^{\top}\Delta y=(x^{\top}\Delta y)^{\top}=y^{\top}\Delta
x=(y,\Delta x)=(\Delta x,y).$

  2. 2.

It is direct to check that

$(\Delta x,x)=\sum_{\\{i,j\\}\in
E}\left(\frac{x_{i}}{\sqrt{d_{i}}}-\frac{x_{j}}{\sqrt{d_{j}}}\right)^{2}\geq
0.$

  3. 3.

It is clear that $\Delta x=0$ if and only if $Ly=0$ where $y=D^{-1/2}x$. Since
the solution of $Ly=0$ is $y=c(1,1,\cdots,1)$ for $c\in\mathbb{R}$, we have
$x=D^{1/2}y=c(\sqrt{d_{1}},\cdots,\sqrt{d_{n}})$, where $c\in\mathbb{R}$.

  4. 4.

The proof is easy and thus we omit it.

∎

The preceding properties have consequences for the eigenvalues of $\Delta$. We
write them as $\lambda_{k}(\Delta)$ so that the eigenvalue equation becomes

$\Delta x^{k}=\lambda_{k}(\Delta)x^{k},$

with $x^{k}$ being a corresponding _eigenvector_.

NB! The eigenvalues and eigenvectors of the normalized Laplacian $\Delta$ can
be rewritten by virtue of the (unnormalized) Laplacian in the following
component-wise form:

$(Lx)_{i}=\lambda\deg i\,x_{i},\;\;i\in V.$

We then order the eigenvalues of $\Delta$ as

$\lambda_{1}(\Delta)=0<\lambda_{2}(\Delta)\leq\cdots\leq\lambda_{n}(\Delta)$

with the $<$ justified by

######  Corollary 12 (Properties of normalized Laplacian matrix spectrum).

  1. 1.

$\lambda_{i}(\Delta)\geq 0,\forall i=1,\ldots,n$.

  2. 2.

The smallest eigenvalue is $\lambda_{1}(\Delta)=0$. Since we assume that $G$
is connected, this eigenvalue is simple, that is $\lambda_{k}(\Delta)>0$ for
$k>1$.

  3. 3.

$\sum_{k=1}^{n}\lambda_{k}(\Delta)=n.$

[1em]0em

###### Proof.

These arguments easily hold from the previous lemma. ∎

######  Lemma 46.

Let $\mathcal{H}^{k}$ be the collection of all $k$-dimensional linear
subspaces of $\mathbb{R}^{n}$. We have the following min-max characterization
of the eigenvalues:

$\lambda_{k}=\max_{H_{k}\in\mathcal{H}^{k-1}}\min\left\\{\frac{(\Delta
x,x)}{(x,x)}:x\neq 0,\ (x,y)=0\mbox{ for all }y\in H_{k}\right\\}$

and dually

$\displaystyle\lambda_{k}$
$\displaystyle=\min_{H_{k}\in{\mathcal{H}^{k}}}\max\left\\{\frac{(\Delta
x,x)}{(x,x)}:x\in H_{k}\backslash\\{0\\}\right\\}$
$\displaystyle=\max_{H_{n-k+1}\in\mathcal{H}^{n-k+1}}\min\left\\{\frac{(\Delta
x,x)}{(x,x)}:x\in H_{n-k+1}\backslash\\{0\\}\right\\}.$

[1em]0em

###### Proof.

This is the well-known min-max theorem for the symmetric matrix $\Delta$. ∎

######  Lemma 47.

  1. (a)

Any eigenvector $x$ for some eigenvalue $\lambda\neq 0$ satisfies

$\lambda=\frac{(\Delta x,x)}{(x,x)}.$

  2. (b)

The second smallest eigenvalue of $\Delta$ is given by

$\lambda_{2}=\min\left\\{\frac{(\Delta x,x)}{(x,x)}:\sum_{i\in
V}\sqrt{d_{i}}x_{i}=0\right\\}=\min\left\\{\left.\frac{\sum\limits_{\\{i,j\\}\in
E}(y_{i}-y_{j})^{2}}{\sum_{i}\deg i\ y_{i}^{2}}\right|\sum_{i}\deg i\
y_{i}=0\right\\}$

  3. (c)

The largest eigenvalue of $\Delta$ is given by

$\lambda_{n}=\max\left\\{\frac{(\Delta
x,x)}{(x,x)}:x\in\mathbb{R}^{n}\setminus\\{0\\}\right\\}.$

  4. (d)

The eigenvalues satisfy

$0\leq\lambda\leq 2.$

[1em]0em

###### Proof.

  1. (a)

It follows from $\Delta x=\lambda x$ that

$\frac{(\Delta x,x)}{(x,x)}=\frac{(\lambda
x,x)}{(x,x)}=\lambda\frac{(x,x)}{(x,x)}=\lambda.$

  2. (b)

Since $(\sqrt{d_{1}},\cdots,\sqrt{d_{n}})$ is an eigenvector corresponding to
$\lambda_{1}=0$, the first equality is due to min-max theorem of Rayleigh
quotient in linear algebra. The second equality is then followed by taking
$x_{i}=\sqrt{d_{i}}y_{i}$.

  3. (c)

This is also due to min-max theorem of Rayleigh quotient.

  4. (d)

Suppose that $x$ is an eigenvector corresponding to $\lambda$, then by (a), we
have

$\lambda=\frac{(\Delta x,x)}{(x,x)}=\frac{\sum_{\\{i,j\\}\in
E}\left(\frac{x_{i}}{\sqrt{d_{i}}}-\frac{x_{j}}{\sqrt{d_{j}}}\right)^{2}}{\sum_{i\in
V}x_{i}^{2}}.$

By Cauchy’s inequality,

$\sum_{\\{i,j\\}\in
E}\left(\frac{x_{i}}{\sqrt{d_{i}}}-\frac{x_{j}}{\sqrt{d_{j}}}\right)^{2}\leq\sum_{\\{i,j\\}\in
E}2\left(\frac{x_{i}^{2}}{d_{i}}+\frac{x_{j}^{2}}{d_{j}}\right)=2\sum_{i\in
V}\sum_{j\in V:j\sim i}\frac{x_{i}^{2}}{d_{i}}=2\sum_{i\in
V}d_{i}\frac{x_{i}^{2}}{d_{i}}=2\sum_{i\in V}x_{i}^{2}$

implying that $\lambda\leq 2$.

∎

######  Theorem 30.

Among all the graphs with $n$ vertices, the complete graph has the largest
possible $\lambda_{2}(\Delta)$ and the smallest possible
$\lambda_{n}(\Delta)$.

[1em]0em

###### Proof.

In fact, since $\lambda_{1}(\Delta)=0$ and the trace
$\sum_{i=1}^{n}\lambda_{i}(\Delta)=n$, we have
$(n-1)\lambda_{2}(\Delta)\leq\sum_{i=2}^{n}\lambda_{i}(\Delta)=n$ which means
$\lambda_{2}(\Delta)\leq n/(n-1)$. The equality holds iff
$\lambda_{2}(\Delta)=\lambda_{3}(\Delta)=\cdots=\lambda_{n}(\Delta)=\frac{n}{n-1}$.
In this case, the eigenvalue $n/(n-1)$ has multiplicity $n-1$, and thus every
vector $x$ satisfying $\sum_{i=1}^{n}x_{i}=1$ must be an eigenvector. In
particular, $(1,-1,0,\cdots,0)$ is an eigenvector corresponding to $n/(n-1)$.
Thus, $1+\frac{1}{d_{1}d_{2}}=(\Delta x,x)/(x,x)=\lambda=\frac{n}{n-1}$, where
we have put $x=(1,-1,0,\cdots,0)$. This derives $d_{1}d_{2}=(n-1)^{2}$, and
then together with $d_{i}\leq n-1$, we have $d_{1}=d_{2}=n-1$. Similarly, we
can obtain $d_{i}=n-1$ for any $i$, which means that $G$ is only possible to
be a complete graph. ∎

NB! This result indicates that complete graphs can be uniquely characterized
by the spectrum of normalized Laplacian. Below, we show a characterization of
multi-partite graphs, meaning that in some sense, multi-partite graphs can be
partially characterized by the spectrum of normalized Laplacian.

######  Definition 16.2.

An undirected graph is called $k$-partite _iff._ one can divide the vertices
by $k$ disjoint groups such that any two vertices from the same group are not
adjacent.

In the same way one can define complete $k$-partite graph:

######  Definition 16.3.

A $k$-partite graph is called complete $k$-partite graph if any two vertices
from the different groups are adjacent.

######  Theorem 31.

A graph with
$\lambda_{t+1}(\Delta)>\lambda_{t}(\Delta)=...=\lambda_{2}(\Delta)=1$ must be
a complete $(n-t+1)$-partite graph, where $2\leq t\leq n-1$.

[1em]0em

###### Proof.

We first prove an elementary observation that for any given $k=1,\cdots,n$,
the following statements hold:

  * •

$\lambda_{k}(\Delta)>1$ if and only if $\lambda_{n-k+1}(A)<0$

  * •

$\lambda_{k}(\Delta)<1$ if and only if $\lambda_{n-k+1}(A)>0$

  * •

$\lambda_{k}(\Delta)=1$ if and only if $\lambda_{n-k+1}(A)=0$

In fact, by the relation $\Delta=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$, it is
clear that
$\lambda_{k}(\Delta)=1-\lambda_{n-k+1}(D^{-\frac{1}{2}}AD^{-\frac{1}{2}})$. By
the min-max theorem of Rayleigh quotient, it is clear that for any
$i=1,\cdots,n$,

$\lambda_{i}(D^{-\frac{1}{2}}AD^{-\frac{1}{2}})=\min_{\dim X\geq i}\max_{x\in
X}\frac{\sum_{\\{u,v\\}\in E}\frac{x_{u}x_{v}}{\sqrt{\deg u\cdot\deg
v}}}{x_{1}^{2}+\cdots+x_{n}^{2}}\mathop{=}\limits^{y_{v}=\frac{x_{v}}{\deg
v}}\min_{\dim X\geq i}\max_{y\in X}\frac{\sum_{\\{u,v\\}\in
E}y_{u}y_{v}}{\sum_{v=1}^{n}\deg v\cdot y_{v}^{2}},$
$\lambda_{i}(A)=\min_{\dim X\geq i}\max_{y\in X}\frac{\sum_{\\{u,v\\}\in
E}y_{u}y_{v}}{\sum_{v=1}^{n}y_{v}^{2}}.$

Since for any $y\in\mathbb{R}^{n}$, the quantity

$\frac{\sum_{\\{u,v\\}\in E}y_{u}y_{v}}{\sum_{v=1}^{n}\deg v\cdot y_{v}^{2}}$

lies between

$\frac{1}{\min_{v}\deg v}\frac{\sum_{\\{u,v\\}\in
E}y_{u}y_{v}}{\sum_{v=1}^{n}y_{v}^{2}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{ and }\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{1}{\max_{v}\deg v}\frac{\sum_{\\{u,v\\}\in
E}y_{u}y_{v}}{\sum_{v=1}^{n}y_{v}^{2}}.$

From this, we can immediately get that $\lambda_{i}(A)>0$ (or $<0$) if and
only if $\lambda_{i}(D^{-\frac{1}{2}}AD^{-\frac{1}{2}})>0$ (or $<0$). Indeed,
all basic properties for the eigenvalue 0 of an adjacency matrix can be
translated into the language of the eigenvalue 1 of the corresponding
normalized Laplacian.

Now, let $G$ be a connected graph which is not complete. If $G$ is not a
complete multipartite graph, then $G$ has three vertices $u,v,w$ such that
$u\sim v$ (i.e., $\\{u,v\\}\in E$), $\\{u,w\\}\not\in E$ and $\\{v,w\\}\not\in
E$. Consider a shortest path from $w$ to $\\{u,v\\}$, we see that $G$ has two
possible types of induced subgraphs:

$\bullet$$w$$\bullet$$\bullet$$\bullet$$\bullet$$u$$v$$w^{\prime}$$w^{\prime\prime}$$\bullet$$w$$\bullet$$\bullet$$\bullet$$\bullet$$u$$v$$w^{\prime}$$w^{\prime\prime}$

Since both of these two graphs on the four vertices
$u,v,w^{\prime},w^{\prime\prime}$, have the second largest eigenvalue larger
than 0, we have $\lambda_{n-1}(A)>0$ by interlacing. Thus,
$\lambda_{2}(\Delta)<1$ which contradicts with the assumption. Hence, $G$ is a
complete multipartite graph.

Suppose that $G$ is a complete multipartite graph on $k$ parts
$V_{1},\cdots,V_{k}$, then we find that $Ax=0$ for any $x$ with the property
that $\sum_{i\in V_{j}}x_{i}=0$, $\forall j=1,\cdots,k$. Thus, 0 is an
eigenvalue and its multiplicity is $n-k$. Let $s_{j}=\sum_{i\in V_{j}}x_{i}$,
$j=1,\cdots,k$. Then $x^{T}Ax=(s_{1}+\cdots+s_{k})^{2}-s_{1}^{2}-\cdots-
s_{k}^{2}>0$ for any positive $x_{i}>0$. Thus, $A$ has exactly one positive
eigenvalue whose multiplicity is 1. Then, we have that
$\lambda_{n-1}(A)=\cdots=\lambda_{n-(n-k)}(A)=0>\lambda_{k-1}(A)$, which
derives that
$\lambda_{n+1-(k-1)}(\Delta)>\lambda_{n+1-k}(\Delta)=\cdots=\lambda_{2}(\Delta)=1$.
Taking $k=n-t+1$, we obtain the desired relation in the theorem. ∎

###  16.2 Cheeger inequality

We consider a subset $S$ of the vertex set $V$, and its complement
$\overline{S}=V\setminus S$. We define the _volume_ of $S$ as

$vol(S)=\sum_{v\in S}d_{v},$

and the _boundary measure_ of $S$ as

$|\partial S|=|E(S,\overline{S})|=\\#\\{\\{i,j\\}\in E;i\in S,j\not\in S\\}.$

The reason that we call the set $\partial S:=\\{\\{i,j\\}\in E;i\in S,j\not\in
S\\}$ the _boundary_ of $S$ is that each edge $\\{i,j\\}$ from $\partial S$
cross both $S$ and $S^{c}$. This property is in the same spirit of boundary in
topology – a point is called a boundary point of a set $S$ if every (closed)
neighbourhood of such point intersects both $S$ and the complement $S^{c}$.

In fact, when we consider a graph $G$ as a one-dimensional simplicial complex,
the smallest closed star neighbourhood of an edge is the collection of its two
end-vertices and such edge itself. So, every neighbourhood of such edge
$\\{i,j\\}\in\partial S$ intersects both $S$ and the complement $S^{c}$.

$\bullet$$S$$S^{c}$$\circ$$\circ$$S$$S^{c}$

NB! Note that in the left picture, in the original topological setting, any
small neighbourhood (marked in orange) of a boundary point intersects the sets
$S$ and $S^{c}$. In the same spirit, on the right picture for the setting of
graph, the smallest closed neighbourhood (marked in orange) of any boundary
edge (marked in orange) has one endpoint in $S$ and the other in $S^{c}$.

We then introduce the _Cheeger constant_

$h=\min_{S\text{ nonempty proper subset of
}V}\frac{|E(S,\overline{S})|}{\min\\{\mathrm{vol}(S),\mathrm{vol}(\overline{S})\\}}$

Now let us state a strong version of the Cheeger inequality:

######  Theorem 32.

$1-\sqrt{1-h^{2}}\leq\lambda_{2}(\Delta)\leq 2h.$

[1em]0em

###### Proof.

Let the vertex set $V$ be divided into the two disjoint sets $U,\overline{U}$
of nodes, and let $U$ be the one with the smaller volume. We consider a vector
$x$ that equals 1 on all the nodes in $U$ and equals $-\alpha$ for some
positive $\alpha$ on $\overline{U}$ so that the equation $\sum_{i\in
V}x_{i}\deg i=0$ holds, that is,

$\sum_{i\in U}\deg i-\alpha\sum_{i\in\overline{U}}\deg i=0.$

Since $\overline{U}$ is the subset with the larger volume
$\sum_{i\in\overline{U}}\deg i$, we have $\alpha\leq 1$. Thus, for our choice
of $x$, the quotient becomes

$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\min\left\\{\left.\frac{\sum\limits_{\\{i,j\\}\in
E}(x_{i}-x_{j})^{2}}{\sum_{i}\deg i\ x_{i}^{2}}\right|\sum_{i}\deg i\
x_{i}=0\right\\}$ $\displaystyle\leq$
$\displaystyle\frac{(1+\alpha)^{2}|E(U,\overline{U})|}{\sum_{i\in U}\deg i\
+\sum_{i\in\overline{U}}\deg i\ \alpha^{2}}$ $\displaystyle=$
$\displaystyle\frac{(1+\alpha)^{2}|E(U,\overline{U})|}{\sum_{i\in U}\deg
i+\sum_{i\in U}\deg i\,\alpha}$ $\displaystyle=$
$\displaystyle\frac{(1+\alpha)|E(U,\overline{U})|}{\sum_{i\in U}\deg i}$
$\displaystyle\leq$ $\displaystyle 2\frac{|E(U,\overline{U})|}{\sum_{i\in
U}\deg i}=2\frac{|E(U,\overline{U})|}{\mbox{vol}(U)}.$

Since this holds for all such splittings of our graph $G$, we obtain the upper
bound

$\lambda_{2}\leq 2h.$

Next, we turn to the hard part. Let $y$ be an eigenvector corresponding to
$\lambda_{2}$. Since $\sum_{i}y_{i}=0$ and $y\neq 0$, we have $\\{i\in
V:y_{i}>0\\}\neq\varnothing$ and $\\{i\in V:y_{i}<0\\}\neq\varnothing$.
Without loss of generality, we may assume that

$\sum_{i\in V:y_{i}>0}d_{i}\leq\frac{1}{2}\sum_{i\in V}d_{i}.$

Let $x$ be the restriction of $y$ onto $\\{i\in V:y_{i}>0\\}$, that is,

$x_{i}=\begin{cases}y_{i},&\text{ for }i\in V\text{ with }y_{i}>0\\\ 0,&\text{
otherwise}.\\\ \end{cases}$

We shall first prove that

$R(x):=\frac{\sum\limits_{\\{i,j\\}\in E}(x_{i}-x_{j})^{2}}{\sum_{i}\deg i\
x_{i}^{2}}\leq\lambda_{2}.$

For $i\in V$ with $x_{i}>0$,

$(Lx)_{i}=\sum_{j\in V:i\sim j}(x_{i}-x_{j})\leq\sum_{j\in V:i\sim
j}(y_{i}-y_{j})=\lambda_{2}d_{i}y_{i}=\lambda_{2}d_{i}x_{i},$

meaning that

$(Lx,x):=\sum_{i\in V}(Lx)_{i}x_{i}=\sum_{i\in
V:x_{i}>0}(Lx)_{i}x_{i}\leq\sum_{i\in
V:x_{i}>0}\lambda_{2}d_{i}x_{i}x_{i}=\lambda_{2}\sum_{i\in V}d_{i}x_{i}^{2}.$

Thus, we have $R(x)\leq\lambda_{2}$.

Also, it follows from

$\frac{\sum\limits_{\\{i,j\\}\in E}(x_{i}+x_{j})^{2}}{\sum_{i}\deg i\
x_{i}^{2}}=2-\frac{\sum\limits_{\\{i,j\\}\in E}(x_{i}-x_{j})^{2}}{\sum_{i}\deg
i\ x_{i}^{2}}=2-R(x)$

that

$\left(\frac{\sum\limits_{\\{i,j\\}\in E}|x_{i}^{2}-x_{j}^{2}|}{\sum_{i}\deg
i\ x_{i}^{2}}\right)^{2}\leq\frac{\sum\limits_{\\{i,j\\}\in
E}(x_{i}-x_{j})^{2}\sum\limits_{\\{i,j\\}\in
E}(x_{i}+x_{j})^{2}}{(\sum_{i}\deg i\ x_{i}^{2})^{2}}=R(x)(2-R(x)).$ (2)

Next we shall use a standard procedure in discrete analysis to show

$\frac{\sum\limits_{\\{i,j\\}\in E}|x_{i}^{2}-x_{j}^{2}|}{\sum_{i}\deg i\
x_{i}^{2}}\geq h.$

In fact, taking $V_{t}=\\{i\in V:x_{i}^{2}>t\\}$ with $t\geq 0$, and taking
$M=\max_{i\in V}x_{i}^{2}$, we have

$\int_{0}^{M}vol(V_{t})dt=\int_{0}^{M}\sum_{i\in
V_{t}}d_{i}dt=\int_{0}^{M}\sum_{i\in V}1_{x_{i}^{2}>t}d_{i}dt=\sum_{i\in
V}\int_{0}^{M}1_{x_{i}^{2}>t}d_{i}dt=\sum_{i\in V}x_{i}^{2}d_{i}$

and

$\displaystyle\int_{0}^{M}|\partial V_{t}|dt$
$\displaystyle=\int_{0}^{M}\sum_{\\{i,j\\}\in E:i\in V_{t},j\not\in
V_{t}}1dt=\int_{0}^{M}\sum_{\\{i,j\\}\in E}1_{x_{i}^{2}>t\geq x_{j}^{2}}dt$
$\displaystyle=\sum_{\\{i,j\\}\in E}\int_{0}^{M}1_{x_{i}^{2}>t\geq
x_{j}^{2}}dt=\sum_{\\{i,j\\}\in E}|x_{i}^{2}-x_{j}^{2}|.$

Now, take $t_{0}\in[0,M]$ such that

$\frac{|\partial V_{t_{0}}|}{vol(V_{t_{0}})}=\min_{0\leq t\leq
M}\frac{|\partial V_{t}|}{vol(V_{t})}.$

Then, $vol(V_{t_{0}})\leq vol(V_{0})=vol(\\{i\in V:x_{i}>0\\})\leq vol(V)/2$
which implies $vol(V_{t_{0}})\leq vol(V\setminus V_{t_{0}})$ and thus

$h\leq\frac{|\partial V_{t_{0}}|}{\min\\{vol(V_{t_{0}}),vol(V\setminus
V_{t_{0}})\\}}=\frac{|\partial V_{t_{0}}|}{vol(V_{t_{0}})}\leq\frac{|\partial
V_{t}|}{vol(V_{t})}$

for any $0\leq t\leq M$. So,

$|\partial V_{t}|\geq vol(V_{t})h$

and hence

$\int_{0}^{M}|\partial V_{t}|dt\geq h\int_{0}^{M}vol(V_{t})dt.$

Therefore, together with all the above facts, we get

$\frac{\sum\limits_{\\{i,j\\}\in E}|x_{i}^{2}-x_{j}^{2}|}{\sum_{i}\deg i\
x_{i}^{2}}=\frac{\int_{0}^{M}|\partial V_{t}|dt}{\int_{0}^{M}vol(V_{t})dt}\geq
h.$

By (2), we have

$R(x)(2-R(x))\geq h^{2}$

which implies

$1-\sqrt{1-h^{2}}\leq R(x)\leq 1+\sqrt{1-h^{2}}.$

We finally obtain

$\lambda_{2}\geq R(x)\geq 1-\sqrt{1-h^{2}}.$

The whole proof is then completed. ∎

Note that $0\leq h\leq 1$, and this yields $0\leq
1-h^{2}\leq(1-\frac{h^{2}}{2})^{2}$ which is equivalent to $h^{2}/2\leq
1-\sqrt{1-h^{2}}$. Then, by the strong Cheeger inequality in Theorem 32, we
immediately obtain the usual Cheeger inequality $h^{2}/2\leq\lambda_{2}$.

##  17 Graph $p$-Laplacian

To illustrate the $p$-Laplacian on graphs, we need to recall the formulation
of the classical $p$-Laplacian $L_{p}$ acting on a smooth function $f$:

$L_{p}f=-\mathrm{div}(|\nabla f|^{p-2}\nabla f)$

where $\nabla$ and $\mathrm{div}$ represent the gradient operator and the
divergence operator, respectively.

Let $G=(V,E)$ be a finite, undirected, simple graph with $n$ vertices and $m$
edges. Along this spirit, we define the gradient $\nabla:C(V)\to C(E)$ as

$(\nabla x)([i,j])=x_{j}-x_{i},$

where

$C(V):=\\{\text{real-valued functions on }V\\}\cong\mathbb{R}^{n}$

and

$C(E):=\\{\text{real-valued functions on }E\\}\cong\mathbb{R}^{m}.$

Now we can think of an undirected graph $G$ as a directed graph with two
orientations, denoted by $[i,j]$ and $[j,i]$, on each edge $\\{i,j\\}$. And we
shall fix an orientation. A vector field $\psi$ is a map $\psi:E\to\mathbb{R}$
with the property $\psi([i,j])=-\psi([j,i])$ for every $i\sim j$. Note that
this concept is similar to the flow introduced in Section 14.

The divergence $\mathrm{div}:C(E)\to C(V)$ is defined by

$\mathrm{div}\psi(i)=\sum_{j\in V:j\sim i}\psi([i,j]).$

We remark here that due to the anti-symmetry of $\psi$, for any nonempty
proper subset $S\subset V$,

$\mathrm{div}\psi(S):=\sum_{\\{i,j\\}\in E:i\in S,j\not\in S}\psi([i,j])$

indicates the ‘rate of flow’ from $S$ to its complement $V\setminus S$. That
is the reason why we call $\mathrm{div}$ the divergence on graphs.

For $p>1$, the graph $p$-Laplacian $L_{p}:C(V)\to C(V)$ is defined by

$L_{p}x=-\mathrm{div}(|\nabla x|^{p-2}\nabla x).$

Clearly, if $p=2$, we obtain $L_{2}x=-\mathrm{div}(\nabla x)$ which is indeed
the standard Laplacian matrix.

Let $\phi_{p}:\mathbb{R}\to\mathbb{R}$ be defined via
$\phi_{p}(t)=|t|^{p-2}t$. Then

$(L_{p}x)_{i}=\sum_{j\in V:j\sim i}\phi_{p}(x_{i}-x_{j}),\;\forall i\in V.$

By setting $\delta=\nabla$ and $\delta^{*}=-\mathrm{div}$, we have

$L_{p}=\delta^{*}\phi_{p}\delta.$

It should be noted that $\delta=B$ is actually the incidence matrix of $G$,
and the adjoint operator $\delta^{*}$ is the transport $B^{\top}$.

###  17.1 Unnormalized version of the $p$-Laplacian eigenproblem

We say $({\lambda},{x})\in\mathbb{R}\times(C(V)\setminus\\{0\\})$ is an
eigenpair, if it satisfies the eigen-equation

$L_{p}{x}={\lambda}\,\phi_{p}(x)$

i.e.,

$\sum_{j\in V:j\sim
i}\phi_{p}(x_{i}-x_{j})={\lambda}\,\phi_{p}(x_{i}),\;\;i\in V.$ (3)

The eigenvalues / eigenfunctions (also called _eigenvectors_) of $L_{p}$ are
the critical values / critical points of the $p$-Rayleigh quotient

$\mathcal{R}_{p}(x)=\frac{\|\nabla
x\|_{p}^{p}}{\|x\|_{p}^{p}}:=\frac{\sum_{\\{i,j\\}\in
E}|x_{i}-x_{j}|^{p}}{\sum_{i\in V}|x_{i}|^{p}}$

For $p=1$, the eigen-equation for $L_{1}$ reduces to the differential
inclusion:

$0\in\partial_{x}\sum_{\\{i,j\\}\in
E}|x_{i}-x_{j}|-\lambda\,\partial_{x}\sum_{i\in V}|x_{i}|.$

And inspired by this inclusion condition, we define

$L_{1}x:=\partial_{x}\sum_{\\{i,j\\}\in E}|x_{i}-x_{j}|=\lim\limits_{p\to
1,\hat{x}\to x}L_{p}\hat{x}$

i.e., the limit points of $L_{p}\hat{x}$ for $p$ tends to 1, and $\hat{x}$
tends to $x$. However, for the sake of brevity in this section, we focus only
on $L_{p}$ when $p>1$.

######  Statement 5.

If $x$ is an eigenvector corresponding to a nonzero eigenvalue of $L_{p}$,
then

$\sum_{i=1}^{n}\phi_{p}(x_{i})=0.$

[1em]0em

###### Proof.

Since $({\lambda},{x})$ is an eigenpair of $L_{p}$ with $\lambda\neq 0$,

$\sum_{j\in V:j\sim
i}\phi_{p}(x_{i}-x_{j})={\lambda}\,\phi_{p}(x_{i}),\;\;i=1,\cdots,n.$

Summing up the above $n$ equations, we get

$\sum_{i\in V}\sum_{j\in V:j\sim i}\phi_{p}(x_{i}-x_{j})={\lambda}\,\sum_{i\in
V}\phi_{p}(x_{i}).$

Note that the left hand side equals

$\sum_{\\{i,j\\}\in
E}\big{(}\phi_{p}(x_{i}-x_{j})+\phi_{p}(x_{j}-x_{i})\big{)}=0$

where we used the easy fact that $\phi_{p}(t)+\phi_{p}(-t)=0$, $\forall
t\in\mathbb{R}$. This implies that the right hand side is zero, and thus by
$\lambda\neq 0$, we have

$\sum_{i\in V}\phi_{p}(x_{i})=0.$

∎

######  Statement 6.

For any eigenvalue $\lambda$ of $L_{p}$, $\lambda\leq 2^{p-1}d$ where $d$ is
the maximum degree of the graph. The equality $\lambda=2^{p-1}d$ holds _iff._
$G$ has a regular bipartite component containing a vertex realizing the
maximum degree.

[1em]0em

###### Proof.

Since $|x_{i}-x_{j}|^{p}\leq 2^{p-1}(|x_{i}|^{p}+|x_{j}|^{p})$, we have

$\sum_{\\{i,j\\}\in
E}|x_{i}-x_{j}|^{p}\mathop{\leq}\limits^{(*)}2^{p-1}\sum_{\\{i,j\\}\in
E}(|x_{i}|^{p}+|x_{j}|^{p})=2^{p-1}\sum_{i\in
V}d_{i}|x_{i}|^{p}\mathop{\leq}\limits^{(**)}d2^{p-1}\sum_{i\in
V}|x_{i}|^{p}.$

The first inequality (*) holds _iff._ $x_{i}=-x_{j}$ for any $\\{i,j\\}\in E$,
_iff._ $G$ has a bipartite connected component.

The second inequality (**) holds _iff._ $d_{i}=d$ for any $i\in V$ with
$x_{i}\neq 0$. Combining the above facts on the equality case for (*) and
(**), we get that the whole equality $\sum_{\\{i,j\\}\in
E}|x_{i}-x_{j}|^{p}=d2^{p-1}\sum_{i\in V}|x_{i}|^{p}$ holds _iff._ $G$ has a
regular bipartite component that contains vertices attaining maximum degree.
It then follows from

$\lambda=\mathcal{R}_{p}(x)=\frac{\sum_{\\{i,j\\}\in
E}|x_{i}-x_{j}|^{p}}{\sum_{i\in V}|x_{i}|^{p}}$

that the desired conclusion in the statement is true. ∎

Finally, we shall compute the $p$-Laplacian eigenvalues of complete graphs.
This was done in a paper by Amghibech [10], but the original computation
contains a small mistake. Here we present a corrected proof.

######  Theorem 33.

Let $G=(V,E)$ be a complete graph with $V=\\{1,\cdots,n\\}$. Then the nonzero
eigenvalues of $L_{p}$ are
$n-\alpha-\beta+(\alpha^{\frac{1}{p-1}}+\beta^{\frac{1}{p-1}})^{p-1}$, where
$\alpha,\beta\in\mathbb{Z}_{+}$ with $\alpha+\beta\leq n$.

######  Statement 7.

Under the above setting, if $A$ and $B$ are two disjoint nonempty subsets of
$V$, then there exist $a>0$ and $b>0$ such that $a\vec{1}_{A}-b\vec{1}_{B}$ is
an eigenvector corresponding to the eigenvalue
$n-|A|-|B|+(|A|^{\frac{1}{p-1}}+|B|^{\frac{1}{p-1}})^{p-1}$ of $L_{p}$.

[1em]0em

###### Proof.

Note that if $a\vec{1}_{A}-b\vec{1}_{B}$ is an eigenvector of $L_{p}$, the
related eigenequation (3) can be written in component-wise form as

$\begin{cases}|B|(a+b)^{p-1}+(n-|A|-|B|)a^{p-1}=\lambda a^{p-1}\\\
|B|b^{p-1}-|A|a^{p-1}=0\\\ |A|(a+b)^{p-1}+(n-|A|-|B|)b^{p-1}=\lambda
b^{p-1}.\end{cases}$

By solving the equations above, we determine $a,b,\lambda$ as follows:
$a=c|B|^{\frac{1}{p-1}}$, $b=c|A|^{\frac{1}{p-1}}$ and
$\lambda=n-|A|-|B|+(|A|^{\frac{1}{p-1}}+|B|^{\frac{1}{p-1}})^{p-1}$, where
$c\neq 0$. ∎

######  Statement 8.

If $x$ is an eigenvector corresponding to a positive eigenvalue of $L_{p}$,
then there exist disjoint nonempty subsets $A,B\subset V$, and $a,b>0$, such
that $x=a\vec{1}_{A}-b\vec{1}_{B}$.

[1em]0em

###### Proof.

Let $(\lambda,x)$ be an eigenpair with $\lambda>0$. Without of loss of
generality, we may assume $x_{1}\leq\cdots\leq x_{n}$ due to the complete
symmetry of $G$.

By Statement 5, there exists $1\leq\alpha<\beta\leq n$ such that

$x_{1}\leq\cdots\leq
x_{\alpha}<0=x_{\alpha+1}=\cdots=x_{\beta-1}=0<x_{\beta}\leq\cdots\leq x_{n}.$

Then, the eigenequation (3) becomes

$\begin{cases}\lambda=-\sum_{k=1}^{i}(-1+\frac{x_{k}}{x_{i}})^{p-1}+\sum_{k=i}^{n}(1-\frac{x_{k}}{x_{i}})^{p-1}&\text{
for }i=1,\cdots,\alpha\\\
\lambda=\sum_{k=1}^{i}(1-\frac{x_{k}}{x_{i}})^{p-1}-\sum_{k=i}^{n}(-1+\frac{x_{k}}{x_{i}})^{p-1}&\text{
for }i=\beta,\cdots,n\\\
\sum_{k=1}^{\alpha}(-x_{k})^{p-1}-\sum_{k=\alpha+1}^{n}x_{k}^{p-1}=0&\end{cases}$
(4)

We shall prove that $x_{1}=\cdots=x_{\alpha}$, and we split the rest of the
proof into two cases:

  * •

$p>2$

In this case, we use the above equations to express $\lambda$ for $i=\alpha$
and for $i=1$.

$\displaystyle\lambda$ $\displaystyle\mathop{=}\limits^{i=1\text{ in
}\eqref{eq:complete:eigen-system}}\sum_{k=1}^{n}(1-\frac{x_{k}}{x_{1}})^{p-1}$
$\displaystyle\leq\sum_{k=1}^{\alpha}(\frac{x_{k}}{x_{\alpha}}-\frac{x_{k}}{x_{1}})^{p-1}+\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{1}})^{p-1}$
(5)
$\displaystyle=\sum_{k=\alpha+1}^{n}(\frac{x_{k}}{x_{1}}-\frac{x_{k}}{x_{\alpha}})^{p-1}+\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{1}})^{p-1}$
(6) $\displaystyle\leq\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{\alpha}})^{p-1}$
(7)
$\displaystyle\leq-\sum_{k=1}^{\alpha}(-1+\frac{x_{k}}{x_{\alpha}})^{p-1}+\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{\alpha}})^{p-1}\mathop{=}\limits^{i=\alpha\text{
in }\eqref{eq:complete:eigen-system}}\lambda$

where the inequality (5) is due to the fact that $1\leq x_{k}/x_{\alpha}$,
$k=1,\cdots,\alpha$, the equality (6) is because of
$\sum_{k=1}^{\alpha}(-x_{k})^{p-1}=\sum_{k=\alpha+1}^{n}x_{k}^{p-1}$, the
inequality (5) is based on the elementary inequality
$a^{t}+b^{t}\leq(a+b)^{t}$ whenever $a,b>0$, $t>1$, in which we take
$a=\frac{x_{k}}{x_{1}}-\frac{x_{k}}{x_{\alpha}}$, $b=1-\frac{x_{k}}{x_{1}}$
and $t=p-1$. So, all the inequalities are in fact equalities, which yield

$x_{1}=\cdots=x_{\alpha}.$

  * •

$1<p<2$

In this case, we use the equations in (4) to represent $\lambda$ when
$i=\alpha$ and $i=n$. We shall use the elementary inequality
$a^{t}\geq(a+b)^{t}-b^{t}$ (or equivalently, $a^{t}+b^{t}\geq(a+b)^{t}$)
whenever $a,b>0$, $0<t<1$, in which we shall take $t=p-1$.

$\displaystyle\lambda$ $\displaystyle\mathop{=}\limits^{i=n\text{ in
}\eqref{eq:complete:eigen-system}}\sum_{k=1}^{n}(1-\frac{x_{k}}{x_{n}})^{p-1}$
$\displaystyle\geq\sum_{k=1}^{\alpha}(\frac{x_{k}}{x_{\alpha}}-\frac{x_{k}}{x_{n}})^{p-1}-\sum_{k=1}^{\alpha}(-1+\frac{x_{k}}{x_{\alpha}})^{p-1}+\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{n}})^{p-1}$
(8)
$\displaystyle=\sum_{k=\alpha+1}^{n}(\frac{x_{k}}{x_{n}}-\frac{x_{k}}{x_{\alpha}})^{p-1}+\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{n}})^{p-1}-\sum_{k=1}^{\alpha}(-1+\frac{x_{k}}{x_{\alpha}})^{p-1}$
(9)
$\displaystyle\geq\sum_{k=\alpha+1}^{n}(1-\frac{x_{k}}{x_{\alpha}})^{p-1}-\sum_{k=1}^{\alpha}(-1+\frac{x_{k}}{x_{\alpha}})^{p-1}\mathop{=}\limits^{i=\alpha\text{
in }\eqref{eq:complete:eigen-system}}\lambda$ (10)

where the inequality (8) is by
$(1-\frac{x_{k}}{x_{n}})^{p-1}\geq(\frac{x_{k}}{x_{\alpha}}-\frac{x_{k}}{x_{n}})^{p-1}-(-1+\frac{x_{k}}{x_{\alpha}})^{p-1}$
for $k=1,\cdots,\alpha$, the equality (9) is because of
$\sum_{k=1}^{\alpha}(-x_{k})^{p-1}=\sum_{k=\alpha+1}^{n}x_{k}^{p-1}$, the
inequality (10) follows from
$(\frac{x_{k}}{x_{n}}-\frac{x_{k}}{x_{\alpha}})^{p-1}+(1-\frac{x_{k}}{x_{n}})^{p-1}\geq(1-\frac{x_{k}}{x_{\alpha}})^{p-1}$.
Again, all the inequalities are actually equalities, which imply
$x_{1}=\cdots=x_{\alpha}.$

If we change $x$ to $-x$, we obtain that $x$ is constant on
$\\{\beta,\beta+1,\cdots,n\\}$; this achieves the proof of the remainder, that
is, we have proved that $x=a\vec{1}_{A}-b\vec{1}_{B}$ for some $a,b>0$, where
$A=\\{1,\cdots,\alpha\\}$ and $B=\\{\beta,\beta+1,\cdots,n\\}$. ∎

Combining Statements 7 and 8, we complete the proof of Theorem 33. It is clear
that when $p=2$, all the positive eigenvalues are equal to $n$.

###  17.2 Normalized $p$-Laplacian eigenvalue problem

The normalized $p$-Laplacian is simply defined by

$\Delta_{p}=D^{-1}L_{p}$

and the relevant eigenproblem is defined as

$\Delta_{p}{x}={\lambda}\,\phi_{p}(x),$

which can be written in the component-wise form like (3) as

$\frac{1}{d_{i}}\sum_{j\in V:j\sim
i}\phi_{p}(x_{i}-x_{j})={\lambda}\,\phi_{p}(x_{i}),\;\;i\in V.$

The spectrum of $\Delta_{p}$ satisfies the following properties.

  1. 1.

If $(\lambda,x)$ is an eigenpair of $\Delta_{p}$, then

$R_{p}(x):=\frac{\sum_{\\{i,j\\}\in E}|x_{i}-x_{j}|^{p}}{\sum_{i\in
V}d_{i}|x_{i}|^{p}}=\lambda$

[1em]0em

###### Proof.

Multiplying $x_{i}$ by

$\sum_{j\in V:j\sim i}\phi_{p}(x_{i}-x_{j})={\lambda}\,d_{i}\phi_{p}(x_{i})$

and summing up the $n$ equalities over $i=1,2,\cdots,n$, we get

$\sum_{i\in V}\sum_{j\in V:j\sim
i}\phi_{p}(x_{i}-x_{j})x_{i}={\lambda}\sum_{i\in V}d_{i}\phi_{p}(x_{i})x_{i}.$

Note that the left hand side is

$\sum_{\\{i,j\\}\in
E}(\phi_{p}(x_{i}-x_{j})x_{i}+\phi_{p}(x_{j}-x_{i})x_{j})=\sum_{\\{i,j\\}\in
E}\phi_{p}(x_{i}-x_{j})(x_{i}-x_{j})=\sum_{\\{i,j\\}\in E}|x_{i}-x_{j}|^{p}$

while the right hand side is

$\lambda\sum_{i\in V}d_{i}\phi_{p}(x_{i})x_{i}=\lambda\sum_{i\in
V}d_{i}|x_{i}|^{p}$

due to the fact that $\phi_{p}(t)t=|t|^{p}$. The proof is then completed. ∎

  2. 2.

The eigenvalues of $\Delta_{p}$ lie in the interval $[0,2^{p-1}]$. [1em]0em

###### Proof.

For any eigenvalue $\lambda$ of $\Delta_{p}$, and for any eigenvector $x$
corresponding to $\lambda$, we have

$0\leq\lambda=\frac{\sum_{\\{i,j\\}\in E}|x_{i}-x_{j}|^{p}}{\sum_{i\in
V}d_{i}|x_{i}|^{p}}\mathop{\leq}\limits^{(*)}\frac{2^{p-1}\sum_{\\{i,j\\}\in
E}(|x_{i}|^{p}+|x_{j}|^{p})}{\sum_{i\in
V}d_{i}|x_{i}|^{p}}=\frac{2^{p-1}\sum_{i\in V}d_{i}|x_{i}|^{p}}{\sum_{i\in
V}d_{i}|x_{i}|^{p}}=2^{p-1}.$

∎

  3. 3.

The largest eigenvalue is $2^{p-1}$ if and only if the graph has a bipartite
connected component. [1em]0em

###### Proof.

The first inequality (*) reduces to an equality _iff._ $x_{i}=-x_{j}$ for any
$\\{i,j\\}\in E$, _iff._ $G$ has a bipartite connected component. ∎

## References

  * [1] Bondy, J. A. (2008). USR Murty Graph Theory. Graduate Texts in Mathematics, 244.
  * [2] XU, Y. (2017). KURATOWSKI’S THEOREM.
  * [3] Patrignani, M. (2013). Planarity Testing and Embedding.
  * [4] Tutte, W. T. (1956). A theorem on planar graphs. Transactions of the American Mathematical Society, 82(1), 99-116.
  * [5] Thomassen, C. (1983). A theorem on paths in planar graphs. Journal of Graph Theory, 7(2), 169-176.
  * [6] KENDALL, M. STEINITZ’THEOREM FOR POLYHEDRA.
  * [7] Strang, A., Haynes, O., Cahill, N. D., Narayan, D. A. (2018). Generalized relationships between characteristic path length, efficiency, clustering coefficients, and density. Social Network Analysis and Mining, 8, 1-6.
  * [8] Chin, C. H., Chen, S. H., Wu, H. H., Ho, C. W., Ko, M. T., Lin, C. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC systems biology, 8(4), 1-7.
  * [9] Cvetković, D., Rowlinson, P., Simić, S. (2009). An Introduction to the Theory of Graph Spectra (London Mathematical Society Student Texts). Cambridge: Cambridge University Press
  * [10] S. Amghibech, Eigenvalues of the discrete $p$-Laplacian for graphs. Ars Comb. 67 (2003), 283-302.