Moore-Penrose Pseudoinverse

Generalization of the inverse of a matrix.

I believe, we pay too much attention to implementation, and too less
attention in the study of the concept that is implemented. I have been
on then teams of many Data Science and Machine Learning projects, and
I would always reiterate on one simple idea; that is, “If you do not
know the math, you don’t know it at all.”

This is a piece of philosophy I deeply believe in. With the advent of packages like numpy, matplotlib, scikit-learn etc., implementing a machine learning model with a moderately difficult data set and problem is fairly simple.
The magic then stays in being able to tweak the algorithm and getting something new (or weird) out of the model. And, for you to be capable of doing so, you will have to know the mechanism behind it.
The Moore-Penrose pseudoinverse in the soul of PCA (Principal Component Analysis), one of the most popularly used Dimensionality reduction techniques.

---

How do we define the inverse of a matrix?
Provided that the matrix is a square matrix and non-singular, we simple divide the adjoint of the matrix with its determinant.
Mathematically, for $A_{m \times m}$ and $|A| \neq 0$, the inverse of $A$ is defined as,
$$A^{-1} = \frac{ \text{adj.} A}{|A|} \tag*{(1)}$$

---

Of course, the above method is computationally very expensive. Hence, we can get the inverse of the matrix recursively using the Fadeev-LeVerrier equation ( Read about that in this blog of mine).

---

Now, how do we deal with matrices that are non-square? How do you find the inverse of a matrix that looks like this,
$$B = \begin{bmatrix} x_{11} & x_{13} \\ x_{21} & x_{23}\\ x_{d1} & x_{d3} \end{bmatrix}_{3\times 2}$$
This is where the Generalization of inverse of a matrix happens, named the Moore-Penrose Pseudoinverse.

---

For every $A_{m \times n}$, there exists a pseudoinverse $A^{\dagger}_{n \times m}$. ( $A^\dagger$ is read as “A dagger”).
$A^{\dagger}$ is mathematically defined as,
$$A^{\dagger}=(A^{T}A)^{-1}A^{T} \tag*{(2)}$$
This is dimensionally consistent. Please check and verify.

---

Now, say we have,
$$Q = \begin{bmatrix} 1 \\ 2 \end{bmatrix}$$
It is impossible to find $Q^{-1}$ by the conventional method $(1)$. So, we use the Generalized Inverse at $(2)$.
So,
$$Q^{T} = \begin{bmatrix} 1 && 2 \end{bmatrix}$$
So, $Q^{T}Q$ comes out to be $[5]_{1 \times 1}$. So, $(Q^{T}Q)^{-1}$ comes out to be $\frac{1}{5}$.
Hence,
$$Q^{\dagger}=(Q^{T}Q)^{-1}Q^{T} =\frac{1}{5} \begin{bmatrix} 1 && 2 \end{bmatrix} = \begin{bmatrix} \frac{1}{5} && \frac{2}{5} \end{bmatrix}$$ which is the pseudoinverse or the generalized inverse.

---

For a square matrix (i.e., $m \times m$),
$$A^{\dagger}=A^{-1}$$
In detail,
$$(A^{T}A)^{-1}A^{T}=\frac{\text{adj. A}}{|A|}$$

---

Some properties of the generalized inverse are,
1. $AA^{\dagger}A=A$
2. $A^{\dagger}AA^{\dagger}=A^{\dagger}$
3. $(AA^{\dagger})=AA^{\dagger}$
One important point to remember is, $A^{\dagger}$ always exists and is unique.

---

Cheers!

The Linear Quadratic Regulator

Optimal Control and Linear-Quadratic-Regulator (LQR)

Today, I will not write an introductory passage to write off my blog. Because, writing an introduction to Optimal Control in itself will required a blog. However, I will add in small tidbits as and when needed.
To understand the topic, we need some basic definitions with us.

---

1.
A control system can be represented in terms of State Space, as follows,
$$\boxed{\dot{x}(t) = A(t)x(t)+B(t)u(t) \\ y(t)=C(t)x(t)+D(t)u(t)}$$
In the above formulation,
$x(.)$ is the state vector; $x(t) \in \mathcal{R}^{n}$.
$y(.)$ is the output vector; $y(t) \in \mathcal{R}^{q}$.
$u(.)$ is the input vector; $u(t) \in \mathcal{R}^{p}$.
$A(.)$ is the System Matrix; $\mathrm{dim}[A(.)]=n\times n$.
$B(.)$ is the Input Matrix; $\mathrm{dim}[B(.)]=n\times p$.
$C(.)$ is the Output Matrix; $\mathrm{dim}[C(.)]=q\times n$.
$A(.)$ is the Feed-forward Matrix; $\mathrm{dim}[D(.)]=q\times p$.
Now, for a system to be controllable, we first define a matrix $\mathcal{Q_c}$, called the controllability matrix, such that,
$$\boxed{\mathcal{Q}_c=\begin{bmatrix} B & AB & A^{2}B & \ldots & A^{n-1}B \\ \end{bmatrix}}$$
The system is controllable if $\mathcal{Q}_c$ has full row rank (i.e. rank($\mathcal{Q}_c$) $=n$).

---

We will assume that we deal with Controllable systems only.

---

Usually, a single input system’s state feedback controller is designed using the Eigen-value method, or Pole Placement method.

---

2.
Pole placement method is the methodology of finding the control vector $\mathcal{U}$ in the form $-\mathcal{k}x$
So, the state space representation changes as,
$$\dot{x}=(A-Bk)x$$
$k$ is found as,
$$|sI-(A-Bk)|=(s-\mu_1)(s-\mu_2)...(s-\mu_n)$$
Here, $\mu$ are the desired pole locations. Note that $k$ is defined as $k=\begin{bmatrix} k_1 & k_2 & k_3 & \ldots & k_n\ \end{bmatrix}$

---

However, for a multi-input system the feedback gain i.e. $k$ is not unique.
Linear Quadratic Control strategy is used to deal with this issue.

---

Now, we dive into the Linear Quadratic Regulator (LQR) formulation, for an $m$-input and $n$-state system with $x \in \mathcal{R}^n$,$u \in \mathcal{R}^m$. Consider a system,
$$\dot{x}=A(t)x(t)+B(t)u(t) \text{ provided }x(0)=x_0 \tag*{(1)}$$
Our aim is to find an open loop control $u(\tau)$, for $\tau \in [t_0, t_f]$ such that we minimize:
$$\boxed{J(u, x_0, t_0, t_f) = \int_{t_0}^{t_f}[x^{T}(t)Q(t)x(t)+u^{T}(t)R(t)u(t)]dt+x(t_f)^{T}Sx(t_f) } \tag*{(2)}$$
where $Q(t)$ and $S$ are symmetric positive semi-definite $n \times n$ matrices.
$R(t)$ is a symmetric positive definite $m \times m$ matrix. Note that $x_0$, $t_0$ and $t_f$ are fixed and given data.
The controller aim is to basically keep $x(t)$ close to 0 especially at $t_f$, which is the final time.
In $(2)$,

• $x^T(t)Q(t)x(t)$ works against the transient response.
• $x^T(t_f)Sx(t_f)$ works against the finite state.
• $u^T(t)R(t)u(t)$ works against the control effort.

The above formulation can regulate the output $y(t) = C(t)x(t)$ near $0$.
Note that, we can define, $S$ and $Q(t)$ as $C^T(t)W(t)C(t)$ where, $W(t) \in R^{r \times r}$
We can now have a theorem as follows,

For a system with fixed initial and final conditions, $\dot{x}=f(x,u,t)$; and clearly $x(t_0)=x_0$. We define our time horizon as $[t_0,t_f]$ such that $t \in [t_0,t_f]$. We find $u(t)$ such that our cost function, $J$ is minimized. $J$ is defined as,
$$J(u(.),x_0)= \phi(x(t_f))+\int_{t_0}^{t_f}L(x(t), u(t), t)dt$$
Here, the first term of $J$ is the final cost and the second term is the recurring cost.

---

Now, we will formulate some important functions that will convert the $J$ which is a constrained optimal control problem to a unconstrained optimal control problem. [THIS MAY NOT MAKE SENSE TO YOU, WHICH IS NATURAL. HOLD ON].
$$\boxed{\dot{\lambda}=-H_x=-\frac{\partial{L}}{\partial{x}}-\lambda^T\frac{\partial{f}}{\partial{x}}} \tag*{(3)}$$
Note that, $\lambda(t)$ ( $\in R^n$ ) is called the Lagrangian.
$H$ is the Hamiltonian operator. Defined in terms of $L$ and $\lambda$ as in $(3)$. Or it can be defined as,
$$\boxed{H(x,u,t) := L(x,u,t)+\lambda^T(t)f(x,u,t)} \tag*{(4)}$$
The above definition is in terms of $L$ as defined in the theorem. So, we define $\lambda$ in the same lines. Just for convenience of computation.
$$\boxed{\lambda^T(t_f)=\frac{\partial \phi}{\partial x}(x(t_f))} \tag*{(5)}$$
$(5)$ can be written as $\lambda^T(t_f)=\phi_x(x(t_f))$
Equation $(3)$,$(4)$ and $(5)$ together form a set of $2n$ differential equations (in $x$ and $\lambda$, obviously) with split boundary conditions at $t_0$ and $t_f$. Now, we can easily define $u(t)$ in terms of $x$ or/and $\lambda$.
As mentioned earlier, the solution is found by converting $J$ from a constrained optimal problem to a constrained optimal problem using a Lagrange multiplier function $\lambda(t)$:
$$\boxed{\bar{J}(u,x_0)=J(u(.),x_0)+\int_{t_0}^{t_f} \lambda^T [f(x,u,t)-\dot{x}]dt} \tag*{(6)}$$
Notice that,
$$\frac{d}{dt}(\lambda^T(t)\dot{x}(t))=\dot{\lambda}^T(t)x(t)+\lambda^T(t)\dot{x} \tag*{(7)}$$
Therefore,
$$\int_{t_0}^{t_f}\lambda^T\dot{x}\ dt=\lambda^T(t_f)\dot{x}(t_f)-\lambda^T(t_0)\dot{x}(t_0)-\int_{t_0}^{t_f}\dot{\lambda}^Tx\ dt \tag*{(8)}$$
As the Hamiltonian Function is defined in $(4)$, thus,
$$\boxed{\bar{J}=\phi(x(t_f))-\lambda^T(t_f)x(t_f) + \lambda^T(t_0)x(t_0)+\int_{t_0}^{t_f}[H(x(t),u(t),t)+\dot{\lambda}(t)x(t)] \ dt} \tag*{(9)}$$
The necessary condition for an optimal solution is $\delta \bar{J}$ of the modified cost with respect to all variations of the system be minimal at all times from $t_0$ to $t_f$.
We will define $\delta \bar{J}$ analytically in the next post and formulate the Riccati Equation that will lay the foundation to some amazing control strategies.
Cheers!

i!

Define the Factorial of a Complex number.

In usual sense, factorial is defined as,
$$n! = \prod_{k=1}^{n} k=1 \cdot2\cdot3\ldots(n-2)\cdot(n-1)\cdot n \tag*{(1)}$$
Now, the not so usual definition is based on the famous Gamma Function,
$$\Gamma(t) = \int_0^\infty e^{-x}x^{t-1}\mathrm{d}x \tag*{(2)}$$
There is an unique and very useful property,
$$\Gamma(n+1)=n! \tag*{(3)}$$
To extend $(3)$ into the complex domain, we will first have to go through Analytic Continuation, please read about Analytic Continuation here.
Therefore, after analytic continuation, we can write it as,
$$z! \overset{\text{def}}{=} \Gamma(z+1) \tag*{(4)}$$
For, $z \in \mathcal C-{-1,-2...}$
So, now, clearly,
$$i!= \Gamma(i+1) \tag*{(5)}$$
By $(3)$
$$\Gamma(i+1)=\int_0^\infty e^{-x} x^{(i+1)-1}\mathrm{d}x \tag*{(6)}$$
Clearly,
$$\Gamma(i+1)=\int_0^\infty e^{-x} x^{i}\mathrm{d}x \tag*{(6)}$$
For easier computation, please catch that,
$$x^i = e^{i \ln x} = \cos ( \ln x)+i \sin( \ln x)$$
Let’s break it down,
$$i! = \int_0^\infty e^{-x} (\cos ( \ln x)+i \sin( \ln x))\mathrm{d}x$$
$$\implies i! = \int_0^\infty e^{-x} \cos ( \ln x)\mathrm{d}x+i\int_0^\infty e^{-x} \sin ( \ln x)\mathrm{d}x$$
If you have reached this far, you obviously know how to solve the above integral.
Cheers!

A new Kilogram? What? How?

A new Kilogram? What? How?

To understand this, let us see how time is defined. We all know the SI unit of time is seconds, and it is defined as the international unit of time, the second, is defined by measuring the electronic transition frequency of caesium atoms. Similarly length (SI unit metre) is defined as the length traveled by light in $\frac{1}{3 \times 10^8}$ seconds or to be more precise $\frac{1}{299792458}$ seconds.
The inspiration to write this blog was derived from Veritasium|How We’re Redefining the kg.
We see that, all these are fundamental quantities are define on a given/fixed standard.
However, the kilogram is set as the weight of a metal cylinder in Paris. That is not that ideal, huh!
So, the NIST is trying to define the unit of mass, i.e., kilogram using some already standardized quantities like the Planck’s constant and the Avogadro’s number. The Planck’s constant takes a dive into advanced applied physics (no, not Quantum physics :p ); therefore, demands my blog’s space and time.
The engineering behind it is pretty simple, to put it in the simplest of the terms;

There is a balancing device, which has a mass unit and a coil
unit. The mass unit is balanced with the magnetic field from
the coil using a motor fixed to it, unit and unless, $F_{magnetic} = F_{gravitational}$. To know more, read Watt Balance or Kibble Balance.

Let’s look into the working now, in this process we will generate some cool equations and in the process will learn some new concepts, as and when required.
First, the watt balance, the principle of operation itself says that,
$$F_{magnetic} = F_{gravitational}$$
Taking into account all the usual scientific nomenclature, we can write,
$$\boxed{m \times g = B \times i \times l} \tag*{(1)}$$
Here, $m$ is the mass, $g$ is the acceleration due to gravity, $B$ is the magnetic flux density (in Tesla($T$)), $i$ is the current in the coil and $l$ length of the conductor in the field.
Equation $(1)$ is for the weighing mode operation of Watt Balance. In which, the weights are matched on both sides.
There is another mode, called the velocity mode in which the mass($m$) is lifted at a height and then the coil is moved back and forth in the magnetic field. This motion induces a voltage, $V$, in the coil. By Faraday’s motional emf expression,
$$\boxed{V = B \times l \times v} \tag*{(2)}$$
Here, $v$ is the velocity of the conductor in the magnetic field.
Now, $(1)$ can be written as,
$$B \times l = \frac{m \times g}{i} \tag*{(3)}$$
and $(2)$ can be written as,
$$B \times l = \frac{V}{v} \tag*{(4)}$$
Equating $(3)$ and $(4)$, we have,
$$\frac{m \times g}{i} = \frac{V}{v}$$
Which can be written as,
$$\boxed{m \times g \times v = V \times i} \tag*{(5)}$$
Interestingly, $m \times g \times v$ is the mechanical power (refer Encyclopedia of Electrochemical Power Sources for more info) and $V \times i$ is the electrical power.
You must all be thinking, how does Planck’s constant come to play. Well, please hold your horses. We are almost there.
To measure $V$ in $(5)$, we go into a concept of superconductivity. Called the Josephson Phenomenon. Please read up on its working here. It is also significant because of being the standard of Voltage.
When DC voltage is applied to a Josephson Junction, the junction experiences an oscillation of frequency(read more),
$$f = \frac{2eV}{h}$$
Here, $f$ is the frequency, $e$ is the elementary charge, and $h$ is the Planck’s constant($6.62607004 \times 10^{-34} m^2 kg \cdot s^{-1}$).
The above equation can be written as,
$$\boxed{V = \frac{h\times f}{2e}} \tag*{(6)}$$
For many Junctions ( say $n$) it is,
$$V = n\frac{h\times f}{2e}$$
The Voltage measure here is accurate to $10^{10}$ parts, refer this.
Now, if we write, $i$ as $\frac{V}{R}$ (where $R$ is the resistance offerered by the junction), then, $(5)$ changes to,
$$\boxed{\frac{h^2\times f^2}{4e^2R} = m \times g \times v} \tag*{(7)}$$
Another question now is how do we measure $R$ ?, for that we will use the idea of Quantum Hall Effect. Quantum Hall effect is the standard for resistance, please refer the paper for more information. But suffice to say, the resistance is defined as,
$$\boxed{R = \frac{1}{k}\frac{h}{e^2}}$$
Here, $k$ is $1, 2, 3, ...$
Please note that, without the integer fraction i.e., $R_k = \frac{h}{e^2}$ is called the von-Klitzing constant, this guy got a Nobel for this. Please read more here.

Using the above equation in $(7)$ we get,
$$\frac{h^2\times f^2}{4e^2 \frac{1}{k}\frac{h}{e^2}} = m \times g \times v$$
Which comes to,
$$h^2\times f^2 = m \times g \times v \times 4e^2 \frac{1}{k}\frac{h}{e^2}$$
$$\implies h\times f^2 \times k = m \times g \times v \times 4$$
$$\implies h = \frac{m \times g \times v \times 4}{f^2 \times k} \tag*{(8)}$$
This is for a single Josephson Junction, for $n$ junctions, it becomes,
$$\implies h = \frac{m \times g \times v \times 4}{f^2 \times k \times n^2}$$
Which can look more elegant if we write it as,
$$\boxed{h = \frac{4}{k \times n^2}\frac{g \times v}{f}m}$$
We have seen that $V$ and $R$ were measured very accurately. Similarly, there is a need that we measure the factors very accurately as well. The $v$ is measured using a Laser Interferometer. The $g$ was measured using a Gravimeter.
So, the scientists in NIST, are just putting in some mass $m$ and get the $h$, and keep tuning the value of $m$ till we get a very accurate $h$.
Cheers!

Study of Brachistochrone Problem.

This blog will deal with one of the most elegant topics in mathematics. The famous Brachistochrone curve. Vsuace did a great job at explaining it in this video. This looks good to me, but it doesn’t feel good; because of course the video does not contain tons of fancy mathematical vocabulary and millions of lines of $\LaTeX$.

To get started, we will just put forward a small fact that, this topic is from a larger field of mathematics called the Calculus of Variations.

Before we get into the Analysis of Brachistochrone Curve (which happens to be the heart of this post), let us brush-up on some concepts.

---

Partial Differentiation (advanced preliminary):

Let us look at how we define the derivative of a given function, say we have $g(x)$ and we have,
$$\boxed{\frac{\mathrm{d}}{\mathrm{d}x}g(x)=g'(x)=\lim_{c \to 0}\frac{f(x+ c)-f(x)}{c}}$$
Provided this limit exists.

Now, for a function of many variables it is not easy to compute the total derivative (the usual derivative). Therefore, we calculate the Partial Derivative. You might have already guessed it by now, we define partial derivative as, say we have a function $h(x,y,z)$, then,
$$\boxed{\frac {\partial}{\partial x}h(x,y,z) = \lim_{c \to 0} \frac{h(x+c, y,z)-f(x,y,z)}{c}}$$
Let us look into a small example now, say we have,
$$f(x,y,z) = x^2y^3+x+y+z+\sin(x) \cos(y)$$
Now,
$$\frac{\partial}{\partial x}f(x,y,z) = 2xy^3+1+0+0+\cos(x) \cos(y)$$
similarly,
$$\frac{\partial}{\partial y}f(x,y,z) = 3x^2y^2+0+1+0-\sin(x) \sin(y)$$
and,
$$\frac{\partial}{\partial z}f(x,y,z) = 0 +0+0+1+0$$

---

Idea of Speed and Time (basic preliminary):

Let us consider, the following scenario, (Made using Geogebra)

We have, in the figure we have denoted the path ACB as $\pi_1$ and the path ADB as $\pi_2$.
Consider that, The time taken for a body to move from A to B using path $\pi_1$ is $t_{AB_{\pi_1}}$ similarly, the time taken for a body to move from A to B using path $\pi_2$ is $t_{AB_{\pi_2}}$.
Also, consider that the distance for $\pi_1$ be $l_1$ and the distance for $\pi_2$ be $l_2$.

Then,
Speed ($\pi_1$) =$v_1$= $\frac{l_1}{t_{AB_{\pi_1}}}$ and Speed ($\pi_2$) =$v_2$ =$\frac{l_2}{t_{AB_{\pi_2}}}$.

3. Idea of Distance(basic preliminary):

Distance between two points $(a,b)$ and $(c,d)$ is given by,
$$l = \sqrt{(a-c)^2+(b-d)^2}$$
In case, the distance is measured from origin ($O(0,0)$) to some point $A(x,y)$, we get,
$$l = \sqrt{(x-0)^2+(y-0)^2} = \sqrt{x^2 + y^2}$$

---

Now that the notion is clear, we shall proceed into the analysis.
[NOTE : Some topics which are new to me (or you) will be explained as and when it is required].

---

Time needed for a body to travel from $A$ to $B$ is given by,
$t_{AB}=\frac{l}{v}$ for a linear path.

For a curve or a non-linear path, we consider piece-wise linear distance ($\mathrm{d}l$) and speed ($v$). We can define, $t_{AB}$ as,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\mathrm{d}l}{v}} \tag*{(1)}$$
Now, we must understand a fact that, all distance is composed of $x$ and $y$ coordinates. We consider the point $A$ as the origin ($O(0,0)$) and $B$ as $(x,y)$. Hence,
$$l = \sqrt{x^2 + y^2}\\ \implies \mathrm{d}l = \sqrt{\mathrm{d}x^2 + \mathrm{d}y^2}\\ \implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})\mathrm{d}x^2}$$
$$\boxed{\implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x} \tag*{(2)}$$

---

Also, since the translation is both in the $x$ and $y$ axes, we can say that the speed is gained by $KE$ and $PE$ equality,
$$KE = PE\\ \implies \frac{1}{2}mv^2 = mgy \\ $$
$$\implies \boxed{v = \sqrt{2gy}} \tag*{(3)}$$

---

Use $(2)$ and $(3)$ in $(1)$, we have,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x}{\sqrt{2gy}}} \tag*{(4)}$$
We write
$$\frac{\mathrm{d}y}{\mathrm{d}x}=y'$$
Now, $(4)$ becomes,
$$t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\mathrm{d}x \tag*{(4.1)}$$
Therefore, the function is,
$$\theta = \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\\ \boxed{\theta = (1 +y'^2 )^{1/2}(2gy)^{-1/2}} \tag*{(5)}$$

---

Please notice that, the integral $t_{AB}$ can be written as,
$$t_{AB} = \int_{A}^{B} f(y, y', x) \mathrm{d}x \tag*{(6)}$$
In our case, $f(y, y', x)$ is $\theta$.

Equation$(6)$ closely follows the The Euler-Lagrange equation.

It is pretty simple to define,

---

The Euler-Lagrange equation

For any $I$ such that,
$$I = \int_{a}^{b}f(y,y',x)\mathrm{d}x$$
Where, $y' = \frac{\mathrm{d}y}{\mathrm{d}x}$, then $I$ has an stationary point, if the Euler-Lagrange Equation, given by,
$$\frac{\partial f}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial f}{\partial y'}=0$$
is satisfied.

Therefore, for our analysis,
$$\boxed{\frac{\partial \theta}{\partial y} - \frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial \theta}{\partial y'}=0} \tag*{(7)}$$

---

Stationary value

This is a value at the stationary point.
A stationary point $a$ is the point at which the first derivative of a function becomes zero,
$$f'(x)|_{x=a} = 0 ; a = \text{Stationary point}$$
Side note :
Find the stationary points of $f(x) = x^3-6x^2+9x$.

---

Beltrami Identity

Our $\theta$ (from $(5)$) is such a cool expression, because $x$ does not appear in that; therefore, of course $\boxed{\frac{\partial \theta}{\partial x} = 0}$; this again leads us to another beautiful form, Beltrami Identity, which is given by,
$$\boxed{\theta - y' \frac{\partial \theta}{\partial y'}=C ; C \text{ is a constant}}.\tag*{(8)}$$

---

We now find, $\frac{\partial \theta}{\partial y'}$,
$$\frac{\partial \theta}{\partial y'}=\frac{\partial}{\partial y'}[(1 +y'^2 )^{1/2}(2gy)^{-1/2}] \\ \implies \frac{\partial \theta}{\partial y'} =(2gy)^{-1/2} \frac{1}{2}(1 +y'^2 )^{-1/2}\cdot 2y'$$
$$\therefore \boxed{\frac{\partial \theta}{\partial y'} = (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y'} \tag*{(9)}$$

---

Use $\theta$ and $\frac{\partial \theta}{\partial y'}$ in equation $(8)$, we get,
$$(1 +y'^2 )^{1/2}(2gy)^{-1/2}-y' \cdot (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y' = C$$
$$\implies \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}-\frac{y'^2}{\sqrt{2gy}\sqrt{1+y'^2}} = C$$
$$\implies \frac{1+y'^2-y'^2}{\sqrt{2gy}\sqrt{1+y'^2}}=C$$
$$\implies \boxed{\frac{1}{\sqrt{2gy}\sqrt{1+y'^2}}=C} \tag*{(10)}$$
Rearranging a little gives,
$$\sqrt{1+y'^2}=\frac{1}{\sqrt{2gy}C}$$
$$\implies \boxed{(1+\frac{\mathrm{d}y}{\mathrm{d}x})y = \frac{1}{2gC^2} = M^2} \tag*{(11)}$$
This is the equation of the cycloid as per [4], I am yet to figure out how this happened.
The solution of $(11)$ can be found using a parametric equation,
$$\boxed{x(t) = \frac{1}{2}M^2(t-\sin t)\\ y(t)=\frac{1}{2}M^2(1-\cos t)} \tag*{(12)}$$
To derive the above equation, please refer Math-stackexchange|Derive the parametric equations of a cycloid.

I have even plotted the solutions at $(12)$, we can see that it is a cycloid,

If you wish to play with the visualisation please go to Cycloid | Parametric @ Desmos by Pragyaditya

Cheers!

---

References:
1. Stationary Points
2. Euler Lagrange Equation
3. Derivation of Beltrami Identity
4. Introduction to calculus of variations
5. Brachistochrone @ Wolfram

A Nice Problem Involving Gamma and Beta Function.

[Playing with Gamma and Beta functions].
Solve :
$$\int_{0}^{\pi/2}(\cos^{2}x-\cos^{3}x) \mathrm{d}x$$

Let,
$$I = \int_{0}^{\pi/2}(\cos^{2}x-\cos^{3}x) \mathrm{d}x\tag*{(1)}$$

Beta function is defined as:
$$\boxed{B(x,y) = \int_{0}^{1}t^{x-1}(1-t)^{(y-1)}\mathrm{d}t \hspace{0.2cm}\forall x, y >0}$$

Now, there are many properties of this beautiful mathematical tool, some of them are,

1. $$B(x,y) = B(y,x)$$ They are symmetric.
2. We can write Beta Functions in terms of Gamma Functions as follows, $$B(x,y) = \frac{\Gamma(x) \cdot \Gamma(y)}{\Gamma(x+y)}.$$
3. We can obviously add factorial into the mix too, $$B(x,y)=\frac{(x-1)!(y-1)!}{(x+y-1)}$$
4. We can define it in terms of trigonometric angles too,
   $$B(x,y)=2 \int_{0}^{\pi/2}(\sin \theta)^{2x-1}(\cos \theta)^{2y-1}\mathrm{d}\theta$$

Now,
$$I = \int_{0}^{\infty}(\cos^{2}x)\mathrm{d}x-\int_{0}^{\infty}(\cos^{3}x) \mathrm{d}x \\ \implies I = I_1 + I_2$$

We will solve both separately,
$$I_1 = \int_{0}^{\infty}(\cos^{2}x)\mathrm{d}x \\ \implies I_1 = \int_{0}^{\infty}(\sin x)^0(\cos x)^2\mathrm{d}x \\ \implies I_1 = \int_{0}^{\infty} (\sin x)^{2.\frac{1}{2}-1}(\cos x)^{2.\frac{3}{2}-1}\mathrm{d}x \\ \implies I_1 = \frac{B(\frac{1}{2},\frac{3}{2})}{2}\\ \implies \boxed{ I_1 = \frac{\Gamma(\frac{1}{2}). \Gamma(\frac{3}{2})}{2\Gamma(\frac{1}{2}+\frac{3}{2})}}$$
Similarly, for the other integral,
$$I_2 = \int_{0}^{\infty}(\cos^{3}x)\mathrm{d}x\\ \implies I_2 = \int_{0}^{\infty}(\sin x)^0 (\cos x)^3 \mathrm{d}x\\ \implies I_2 = \int_{0}^{\infty}(\sin x)^{2.\frac{1}{2}-1}(\cos x)^{2.2-1}\mathrm{d}x\\ \implies I_2 = \frac{B(\frac{1}{2},2)}{2}\\ \implies \boxed{ I_2 = \frac{\Gamma(\frac{1}{2}). \Gamma(2)}{2\Gamma(\frac{1}{2}+2)}}$$
Therefore,
$$\boxed{\int_{0}^{\infty}(\cos^{2}x-\cos^{3}x) \mathrm{d}x = \frac{\Gamma(\frac{1}{2}). \Gamma(\frac{3}{2})}{2\Gamma(\frac{1}{2}+\frac{3}{2})}+\frac{\Gamma(\frac{1}{2}). \Gamma(2)}{2\Gamma(\frac{1}{2}+2)}}$$
Cheers $$\bowtie$$

Integration : Problem 2

$$\text{I have been inactive for quite a long time.}\\ \text{As usual, I will start this blogging streak}\\ \text{with a small warm-up exercise.}\\ \text{I will solve some mildly complex Integrals, }$$

$\underline{\text{Question 1.}} \\$ $\LARGE \int _{0}^{\frac{1}{2}}\frac{x^2}{(1+x^2)^{3/2}} \mathrm{d}x$
$\underline{\text{Solution}}$:
$$\text{Let, } I_{1} = \int _{0}^{\frac{1}{2}}\frac{x^2}{(1+x^2)^{3/2}} \mathrm{d}x \\ \text{We can write, } I_1 = \int _{0}^{\frac{1}{2}} \frac{x^2 + 1 - 1}{(1+x^2)^{3/2}} \mathrm{d}x \\ \implies I_1 = \int _{0}^{\frac{1}{2}} \frac{x^2+1}{(1+x^2)^{3/2}}\mathrm{d}x - \int_{0}^{\frac{1}{2}}\frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Which with some algebraic manipulation comes to, }\\ I_1 = \int_{0}^{\frac{1}{2}} \frac{1}{\sqrt{1+x^2}}\mathrm{d}x - \int_{0}^{\frac{1}{2}}\frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Let,} I_2 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x \\ \text{and, Let , } I_3 = \int \frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{We know that, } \boxed{\boxed{\int \frac{\mathrm{d}x}{\sqrt{x^2 + a^2}} = \log|x +\sqrt{x^2 + a^2}|+c }} \\ \therefore I_2 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x = \log|x+\sqrt{1+x^2}| \\ \text{ neglecting } c \text{ because it will ultimately} \\ \text{ be a definite integral }.\\ \\ \text{Now the rigor starts, } I_3 = \int \frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Take, } x = \tan {t} \rightarrow \frac{\mathrm{d}x}{\mathrm{d}t} = \sec^2 {t} \\ \therefore I_3 = \int \frac{\sec^2{t}}{(1+\tan^2{t})^{3/2}} \mathrm{d}t\\ \text{We know, } \boxed{\boxed{1 + \tan^2 \theta = \sec^2 \theta}} \\ \text{Using this } I_3 \text{ changes to, }\\ I_3 = \int \frac{1}{\sec t} \mathrm{d} t \implies \int \cos {t} \mathrm{d}t = \sin t \\ \text{We get the substitution in place ,}\\ t = \tan^{-1} x \rightarrow \sin t \rightarrow \sin{(\tan^{-1} x)} = \frac{x}{\sqrt{x^2+1}} \\ \therefore I_3 = \sin t = \frac{x}{\sqrt{x^2+1}} \text{Now,} \\ I_1 = I_2 - I_3 \text{ with no limits. }\\ \therefore I_1 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x = \log|x+\sqrt{1+x^2}| - \frac{x}{\sqrt{x^2+1}}\\ \text{Put limits } \frac{1}{2} \text{ to } 0 \text{ and get answer as, }\\ I_1 = {I_1}|_{0}^{1/2}=\boxed{0.033}.\\ \text{Cheers! } \bowtie$$