Least-Squares and Solution Uniqueness

I recently was in discussion and could not answer off the top for what to me sounded like a simple question on linear algebra and least squares fitting. So I decided to do some reading and write up a short post on the topic.The question was: is the fitted parameter vector guaranteed to be unique?

For ordinary least squares, the answer is not always. If the design matrix $\mathbf{A}$ has linearly dependent columns, then different parameter vectors $\mathbf{c}$ can produce the same predictions. A common fix is ridge regularization. When the ridge is added with a positive value, the answer becomes yes. The addition of regularization makes the normal-equation matrix positive definite. Positive definite matrices imply invertibility of that linear system, and consequently make the least-squares objective strictly convex, so a unique global minimizer exists.

Just restating it so I remember, when you add a constant diagonal term $\lambda\mathbf{I}$ into the normal-equation matrix (see eq.\ref{eq:H}), you make it positive definite and therefore invertible, which means uniqueness.

The Regularized Problem

Consider the regularized least-squares objective below where we want to find the set of optimal parameters $\mathbf{c}$:

$$ \begin{align} J(\mathbf{c}) &= \left| \mathbf{W}^{1/2} (\mathbf{A}\mathbf{c} - \mathbf{b}) \right|_{2}^{2} \nonumber\\ &\quad + \underbrace{\lambda \left| \mathbf{c} \right|_{2}^{2}}_{\text{ridge regularization term}} \nonumber\\ &\quad + \underbrace{\gamma \left| \mathbf{D}\mathbf{c} \right|_{2}^{2}}_{\text{roughness regularization term}}, \label{eq:objective} \end{align} $$

where $\mathbf{A}$ is the design matrix, $\mathbf{c}$ contains the fitted coefficients, $\mathbf{W}$ is a symmetric positive-semidefinite weight matrix (typically diagonal with nonnegative entries), and $\mathbf{D}$ is a discrete derivative / roughness operator (first- or second-difference, depending on whether one penalizes slope or curvature). Taking the gradient of eq.\ref{eq:objective} and setting it equal to zero yields the normal equation:

$$ \begin{equation} \mathbf{H}\mathbf{c}=\mathbf{A}^{T}\mathbf{W}\mathbf{b}, \label{eq:normal} \end{equation} $$

with $\mathbf{H}$ being defined as:

$$ \begin{equation} \mathbf{H}=\mathbf{A}^{T}\mathbf{W}\mathbf{A}+\lambda\mathbf{I}+\gamma\mathbf{D}^{T}\mathbf{D}. \label{eq:H} \end{equation} $$

The whole question of uniqueness therefore reduces to whether $\mathbf{H}$ is positive definite.

Why is $\mathbf{H}$ Positive Definite?

Let's try to understand this geometrically. Take any nonzero vector $\mathbf{z} = (z_1, z_2, \ldots)$. If we want to see whether a matrix like $\mathbf{H}$ is positive definite, we look at its so-called quadratic form $\mathbf{z}^{T}\mathbf{H}\mathbf{z}$, with $\mathbf{H}$ from eq.\ref{eq:H}. Expanding that definition gives

$$ \begin{equation} \mathbf{z}^{T}\mathbf{H}\mathbf{z}=\left|\mathbf{W}^{1/2}\mathbf{A}\mathbf{z}\right|_{2}^{2}+\lambda\left|\mathbf{z}\right|_{2}^{2}+\gamma\left|\mathbf{D}\mathbf{z}\right|_{2}^{2}. \label{eq:quad} \end{equation} $$

This just means we plug $\mathbf{z}$ into eq.\ref{eq:H} (so, multiply $\mathbf{z}$ by $\mathbf{H}$, then by $\mathbf{z}$ again) and rewrite the result as a sum of squared norms. If the right-hand side of eq.\ref{eq:quad} is always greater than zero for any nonzero vector $\mathbf{z}$, then $\mathbf{H}$ is positive definite.

Figure 1. Equal-cost valleys without ridge (vertical lines; many (z1,z2) pairs share one cost) become closed equal-cost contours with ridge (cost grows even along z1=0).

To try and geometrically visualize eq.\ref{eq:quad}, Figure 1 animates the components of a coefficient vector $\mathbf{z}=(z_{1},z_{2})$. Without ridge, equal-cost contours of the toy form $\mathbf{z}^{T}\mathbf{H}_{0}\mathbf{z}=z_{1}^{2}$ are vertical lines. This means if you fix $z_{1}$ and slide $z_{2}$, the cost does not change. Many different pairs $(z_{1},z_{2})$ therefore share the same cost (i.e., are degenerate solutions). The worst valley is $z_{1}=0$, where $\mathbf{A}\mathbf{z}=\mathbf{0}$ and the cost is zero for every $z_{2}$. When we use regularization (i.e., a ridge), the cost becomes $z_{1}^{2}+\lambda z_{2}^{2}$ and the contours close into a bowl. Even along $z_{1}=0$ one pays $\lambda z_{2}^{2}>0$ unless $\mathbf{z}=\mathbf{0}$ ¹. In jargon, the line $z_{1}=0$ is the null space of $\mathbf{A}$. Therefore,

$$ \mathbf{z}^{T}\mathbf{H}\mathbf{z}>0\qquad\text{for every }\mathbf{z}\neq\mathbf{0}, $$

which means

$$ \mathbf{H}\underbrace{\succ}_{\text{is positive definite}}0. $$

This immediately means that $\mathbf{H}$ is invertible and therefore eq.\ref{eq:normal} has exactly one solution (a unique critical point). Another note is that the Hessian of eq.\ref{eq:objective} is $2\mathbf{H}$, which is also positive definite, so $J$ is strictly convex and that critical point is the unique global minimizer [1].

Where the guarantee actually comes from

The roughness term $\mathbf{D}^{T}\mathbf{D}$ is generally positive semidefinite, not positive definite, because a discrete derivative operator normally has a null space (for example constant vectors for a first-difference $\mathbf{D}$). Roughness regularization alone does not automatically guarantee uniqueness. The unconditional guarantee here comes from the positive ridge term $\lambda\mathbf{I}$; the roughness term adds smoothness without changing that conclusion. Data weights and $\gamma$ do not carry the uniqueness guarantee either.

Closing things out, the compact mathematical statement is then:

$$ \begin{align} \lambda>0 &\Longrightarrow\mathbf{H}\succ0 \nonumber\\ &\Longrightarrow J(\mathbf{c})\text{ is strictly convex} \nonumber\\ &\Longrightarrow\mathbf{c}^{*}\text{ is the unique global minimizer}. \end{align} $$

Punchline

Invertibility of $\mathbf{H}$ is what makes the linear system have one solution. Positive definiteness is what makes $J$ strictly convex, so that solution is the unique global minimizer ... the existence of $\mathbf{H}^{-1}$ is a consequence, not the cause. Numerically you would probably still solve eq.\ref{eq:normal} as a linear system (e.g., via Cholesky decomposition) rather than explicitly forming $\mathbf{H}^{-1}$.

Footnotes

---

1. Same cost is not the same thing as same length $|\mathbf{z}|_{2}$. ↩

References

[1] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. https://doi.org/10.1017/CBO9780511804441.

Reuse and Attribution

Dual Numbers

A while back when I was doing some exploration of writing a simple NN code to improve my understanding of neural networks and deep learning in general, I came across dual numbers. They're a type of number that generalizes the concept of real and complex numbers. But what makes them so interesting is that they can encode both a function value and its derivative in a single number. This means that we can use them to simplify the calculation of derivatives and solve complex problems efficiently.

So how does one think of dual numbers? What's the difference between a dual number and a complex number? One way to think about dual numbers is that they consist of two parts: a scalar part and a skew part. The scalar part is just a regular real number, while the skew part is a multiple of a new number, often denoted as $\epsilon$, that satisfies the property $\epsilon^2=0$. This means that every dual number can be written as $a+b\epsilon$, where $a$ and $b$ are real numbers.

What's most interesting is that the skew part of a dual number is that it provides an approximation of the first derivative of a function evaluated at a particular point. By using the dual number representation of the function at that point, one can calculate both the function value and its derivative in one shot.  One reason dual numbers have applications in deep learning is that algebra on dual numbers provides the chain rule for calculus, therefore they can be used to compute derivatives of complicated functions involving multiple variables and interdependencies.

As an example, say I want to evaluate the function $f(x)=x^2+2x$ at $x=3$. The dual number representation of $f(3)$ is $f(3+\epsilon)=f(3)+f'(3)\epsilon$, where $f'(x)=\frac{df(x)}{dx}$. We can compute $f(3)$ directly as $f(3)=3^2+2\cdot3=9+6=15$. To compute $f'(3)$, we can take the derivative of $f$ with respect to $x$: $f'(x)=2x+2$. Evaluating this at $x=3$, we get $f'(3)=2\cdot3+2=8$. Therefore, the dual number representation of $f$ at $x=3$ is $15+8\epsilon$.

One of the benefits of dual numbers is the derivative of the composition of two functions, $f(g(x))$ requires only the derivatives of the individual functions. Specifically, if $f(x)$ and $g(x)$ are two functions, then the dual number representation of their composition $f(g(x))$ is $(f(g(x)), f'(g(x))g'(x))$. This is especially useful when dealing with complex functions involving multiple variables and complicated interdependencies.

Dual numbers are actually a useful mathematical concept because they have practical applications in a wide range of fields. It's pretty cool that one can encode function values and derivatives in a single number, which makes it possible to simplify the calculation of derivatives and solve complex problems efficiently. On my computational blog, I have an example using dual numbers to calculate the derivative of an interatomic potential, Dual Numbers Pluto blog.

Reuse and Attribution

Simple Example of Entangled States & My Thoughts

Entanglement

In the famous EPR paper, the primary argument is related to the bizarre non-local consequences of entangled quantum states. Einstein's main objection is related to the non-local (spatial) characteristic of entangled states. He suggested that quantum mechanics must be an incomplete theory and argued that our experience in nature appears local. He therefore concluded that hidden variables or additional degrees-of-freedom most likely exist and are not captured by current quantum theory formulations. Unfortunately, well maybe fortunately, Einstein appears to be incorrect about this as exemplified by the work of John S. Bell. The main outcome of Bell's work is that non-local theories can exist without the need for hidden variables and thus entangled quantum states are not restricted by spatial locality. I won't dive to deep into the foundations of quantum theory, but rather just look at a simple mathematical argument for why entangled states exist.

Simple Example

Let me first start with describing a quantum state (i.e. wavefunction) with basis vectors$^*$:

$$|0\rangle = \begin{bmatrix}
1 \\
0
\end{bmatrix} ,\;
|1\rangle = \begin{bmatrix}
0 \\
1
\end{bmatrix}
$$

This is a orthogonal basis set and corresponds to a Hilbert space (abstract function vector space) of $2^n$, where $n$ is the number of particles or objects. Each particle or object can be represented by a linear combination vectors in the Hilbert space. The basis corresponds to that used to describe quantum bits or qubits. 

The basis states can be used to form product basis states, which are given by the tensor product, such that for a system of two objects the  Hilbert space  is $2^2 = 4$ and can be written as:

$$
|00\rangle =\begin{bmatrix}
1 \\
0 \\
0 \\
0 \end{bmatrix} ,\;
|01\rangle =\begin{bmatrix}
0 \\
1 \\
0 \\
0 \end{bmatrix} ,\;
|10\rangle =\begin{bmatrix}
0 \\
0 \\
1 \\
0 \end{bmatrix} ,\;
|11\rangle =\begin{bmatrix}
0 \\
0 \\
0 \\
1 \end{bmatrix} \;,
$$

where the state of a particle/object is represented by a linear superposition of the product basis states, each having a complex amplitude (can be real if the imaginary part is zero). Now we can do something interesting, what if we take the product basis just written above and construct a  potential wavefunction with the following features:

$$ |\Psi\rangle = \frac{1}{\sqrt{2}} \left( |00\rangle + |11\rangle\right).$$

Since we know from that the states $|00\rangle$ and $|11\rangle$ can also be written as product states, as shown above, lets do the following:

$$|\psi_1\rangle = \alpha_1 |0\rangle + \beta_1|1\rangle $$
$$|\psi_2\rangle = \alpha_2 |0\rangle + \beta_2|1\rangle. $$

Now from these definitions lets take the tensor product:

\begin{align}
|\psi\rangle &= |\psi_1\rangle \otimes |\psi_2\rangle \\
 &= \alpha_{1} \alpha_{2} |0\rangle|0\rangle + \alpha_1 \beta_2 |0\rangle|1\rangle + \beta_1 \alpha_2 |1\rangle|0\rangle + \beta_1 \beta_2 |1\rangle|1\rangle \\
 &= \alpha_{1} \alpha_{2} |00\rangle + \alpha_1 \beta_2 |01\rangle + \beta_1 \alpha_2 |10\rangle + \beta_1 \beta_2 |11\rangle .\\
\end{align}

So we now have the given state $|\Psi\rangle$ and the product state $|\psi\rangle$, and if we compare the terms we immediately observe that:

$$\alpha_1 \alpha_2 = \frac{1}{\sqrt{2}} \; \text{and} \;  \beta_1 \beta_2 = \frac{1}{\sqrt{2}}$$

However, since for an orthonormal basis we must have that, $||\Psi\rangle|^2 = 1$ , it then has to be such that:

$$\alpha_1 \beta_2 = 0 \; \text{and} \; \beta_1 \alpha_2 = 0.$$

This would be a contradiction though because its is not feasible  given $\alpha_1 \alpha_2 = \frac{1}{\sqrt{2}}$ and $\beta_1 \beta_2 = \frac{1}{\sqrt{2}}$. We therefore say that quantum state,

$$ |\Psi\rangle = \frac{1}{\sqrt{2}} \left( |00\rangle + |11\rangle\right),$$

represents a quantum entangled state, and cannot be written in terms of product states. Moreover, the determination or expectation of a single particle/object in this entangled state immediately tells us about the state of the other particle/object. For example, if the expectation value of particle 1 is $0.5$ in $|0\rangle$ or $|1\rangle$ then we will know with unity the state of particle 2 is. 

What does it mean

The feature of entanglement is exclusive to quantum mechanics. One way I try to think of things is that quantum states can be composed of two types, those that exist due to the "combination" of quantum states (i.e., product states) and states which manifest as unique quantum solutions which are indecomposable into any other state or description. This fact, of entanglement is one of the great mysteries in quantum physics but its existence is enabling fantastic technologies.

There are two quotes from Niels Bohr's that I think fit well with this post, the first is:

"Einstein, stop telling God what to do [with his dice]!"

-Niels Bohr's, A response to Einstein's assertion that "God doesn't play dice".

This was in response to Einstein's dislike for many of the unexplained conundrums of quantum mechanics. The second quote:

"Anyone who is not shocked by quantum theory has not understood it."

-Niels Bohr's, The Philosophical Writings of Niels Bohr (1987),

Which captures the unexpected and possibly bizarre way of thinking one needs to succumb to in order to appreciate the predictive power of quantum theory (i.e. mathematics and interpretations). Despite this lack of comfort, quantum theories have made extremely accurate predictions and have been validated numerous times through meticulously controlled experiments.

I personally still find the foundations of quantum physics to be nebulous, but this is probably due to my own fallibility. To me its is unsatisfying that we do not not know the true meaning of the wavefunction or more specifically what is the meaning of a Universal wavefunction? Maybe its to complex that we will never know. Then there are questions about why entanglement exist, is it necessary to be consistent with physics as a whole (i.e., General Relativity)? Have we dismissed other understandings to quickly? Does the dendritic many-worlds interpretation originally proposed by Hugh Everett describe reality? What about revisiting non-local hidden variable theories such as Bohmian mechanics (also known as pilot-wave theory and De Broglie-Bohm theory)? How does non-locality make sense, is our notion of space being innate to the universe incorrect? Does space emerge from something else?

All of these question intrigue me, however, I have only scratched the surface  and look forward to learning more about research focused on the foundations of quantum physics. There is a considerable learning curve and start-up time  for me since my formal training is not in theoretical physics, but that won't stop me.

$^*$ The basis vectors are represented using the bra―ket notation pioneered by Paul Dirac. This is a very useful notation but may be unfamiliar to materials science people given that our solid state physics education, to my knowledge, never goes over this because we typically always deal with quantum states (i.e., wavefunctions) in a position or wave-vector basis, for example, $\psi_i(x) = A_i e^{\alpha_i \left(x-x_o\right)^2}$ or $\psi_{k}(x) = \frac{1}{\sqrt{V}} e^{ik\cdot x}$. Therefore the more broad and useful bra―ket notation goes unused.

Reuse and Attribution

Matrices

A matrix is an array of numbers arranged with $m$-rows and $n$-columns. When $m=n$ the matrix is termed square. The index of an element in a matrix is notated using the indices $i$ and $j$ for the rows and columns, respectively. For example:

$$\mathbf{B}=
\begin{bmatrix}
b_{11} & b_{12} \\
b_{21} & b_{22} \\
\end{bmatrix}.
$$

The algebra of matrices is straightforward, for example, addition occurs by adding same indexed elements:

$$\mathbf{B}+\mathbf{C} =
\begin{bmatrix}
b_{11}+c_{11} & b_{12}+c_{12} \\
b_{21}+c_{21} & b_{22}+c_{22} \\
\end{bmatrix}
.$$

We can write the addition operation in more compact form using implicit index notation, e.g., $\mathbf{B}+\mathbf{C}=b_{ij}+c_{ij}$. Keep in mind that the dimensions of the matrices must be the same for addition or subtraction operations.

Multiplication of matrices by a scalar quantity is simply, $S\cdot\mathbf{B}$. When multiplying two matrices, the inner dimensions must be the same, for example, if $\mathbf{B}$ is a $m \times n$ then $\mathbf{C}$ must be a $n \times p$. In other words, the number of columns in $\mathbf{B}$ must be equal to the number of rows in $\mathbf{C}$. The operation is written in compact as follows:

$$ c_{jk} = \sum_{i=1}^{n} a_{ji}b_{ik}.$$

Here are some examples, $c_{11} = a_{11}b_{11}+a_{12}b_{21}$ and $c_{12} = a_{11}b_{12} + a_{12}b_{22}$. An important operation/transformation of matrices is the transpose, which is the process of switching elements in the rows and columns. The transpose is indicated with a superscript capital "T", e.g., $\mathbf{B}^T$. For a square matrix, the diagonal components are commonly referred to as the principal terms and the sum of them is the trace of the matrix. For square matrices that have ones as diagonal terms and zeros as off-diagonal terms, we refer to them as unit of identify matrices, for example:

$$\mathbf{I}=\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{bmatrix}.$$

Identity matrices are commonly represented with $\mathbf{I}$. Another important property of square matrices is invertibility. A matrix is said to have an inverse if it satisfies the following condition:

$$\mathbf{A}\mathbf{B} = \mathbf{I},$$

we call $\mathbf{B}$ the inverse of $\mathbf{A}$. If a square matrix doesn't have an inverse it is referred to as singular. Matrices that have the inverse which satisfies the condition:

$$\mathbf{A}^{T}\mathbf{A}=\mathbf{I},$$

are orthogonal matrices. A determinant is a specific value that can be computed for square matrices and can be thought of as the scaling factor for linear transformations. The determinant is typically written as:

$$\text{det}\, \mathbf{B} = \Delta \mathbf{B} = |\mathbf{B}|,$$

For 2x2 matrices the determinant is determined by taking the product of the diagonals and subsequent difference between them, as shown below:

$$|\mathbf{B}| =
\begin{vmatrix}
b_{11} & b_{12} \\
b_{21} & b_{22} \\
\end{vmatrix} = b_{11}b_{22}-b_{21}b_{12}.$$

Another method for finding the determinant is using the Laplace expansion method. Determinants have the following commuting properties:

$$|\mathbf{A}\mathbf{B}| = |\mathbf{A}||\mathbf{B}|,$$
$$|\mathbf{A}| = |\mathbf{A}^{T}|.$$

A matrix that is singular will have a determinant that evaluates to zero, thus this is a good way for identifying if a matrix has an inverse.

Matrices and determinants are an important component to the field of  linear algebra, so having good command over them is a serious advantage. I know there are a lot of keywords presented in this post that were not given a formal mathematical description but I strongly suggest clicking the link and reading on to get a better sense of there use in linear algebra.

For this post we will provide a quote from a mathematician who did a lot of research in advance algebra and group theories.

"We [Kaplansky and Halmos] share a philosophy about linear algebra: we think basis-free, we write basis-free, but when the chips are down we close the office door and compute with matrices like fury."
-Irving Kaplansky, Paul Halmos: Celebrating 50 Years of Mathematics

References & Additional Reading

[1] Advanced Engineering Mathematics, Wiley, 8th ed., E. Kreyszig.

[2] The Mathematics Companion, IoP, 1st ed., A.C. Fischer-Cripps.

Reuse and Attribution