Logistic Regression using R - How to handle class bias?

Doing this blog as a revision into some primitive topics in Machine Learning.

Since, almost all my current work is done using R, I will implement logistic regression using this language.

We will be using adult data set. After all is set, now we dive straight into the code. Further explanations shall not be provided as most of the code is self-explanatory in nature.

> idata <- read.csv("adult_csv.csv")

> summary(idata)

      age                   workclass         fnlwgt               education     education.num  
 Min.   :0.000   Private         :33906   Min.   :  12285   HS-grad     :15784   Min.   : 1.00  
 1st Qu.:1.000   Self-emp-not-inc: 3862   1st Qu.: 117551   Some-college:10878   1st Qu.: 9.00  
 Median :2.000   Local-gov       : 3136   Median : 178145   Bachelors   : 8025   Median :10.00  
 Mean   :1.771                   : 2799   Mean   : 189664   Masters     : 2657   Mean   :10.08  
 3rd Qu.:3.000   State-gov       : 1981   3rd Qu.: 237642   Assoc-voc   : 2061   3rd Qu.:12.00  
 Max.   :4.000   Self-emp-inc    : 1695   Max.   :1490400   11th        : 1812   Max.   :16.00  
                 (Other)         : 1463                     (Other)     : 7625                  
               marital.status            occupation            relationship                   race      
 Divorced             : 6633   Prof-specialty : 6172   Husband       :19716   Amer-Indian-Eskimo:  470  
 Married-AF-spouse    :   37   Craft-repair   : 6112   Not-in-family :12583   Asian-Pac-Islander: 1519  
 Married-civ-spouse   :22379   Exec-managerial: 6086   Other-relative: 1506   Black             : 4685  
 Married-spouse-absent:  628   Adm-clerical   : 5611   Own-child     : 7581   Other             :  406  
 Never-married        :16117   Sales          : 5504   Unmarried     : 5125   White             :41762  
 Separated            : 1530   Other-service  : 4923   Wife          : 2331                             
 Widowed              : 1518   (Other)        :14434                                                    
     sex         capitalgain      capitalloss      hoursperweek         native.country    class      
 Female:16192   Min.   :0.0000   Min.   :0.0000   Min.   :0.000   United-States:43832   <=50K:37155  
 Male  :32650   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:2.000   Mexico       :  951   >50K :11687  
                Median :0.0000   Median :0.0000   Median :2.000                :  857                
                Mean   :0.2003   Mean   :0.1149   Mean   :1.951   Philippines  :  295                
                3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:2.000   Germany      :  206                
                Max.   :4.0000   Max.   :4.0000   Max.   :4.000   Puerto-Rico  :  184                
                                                                  (Other)      : 2517  

> head(idata)

  age        workclass fnlwgt education education.num     marital.status        occupation  relationship  race
1   2        State-gov  77516 Bachelors            13      Never-married      Adm-clerical Not-in-family White
2   3 Self-emp-not-inc  83311 Bachelors            13 Married-civ-spouse   Exec-managerial       Husband White
3   2          Private 215646   HS-grad             9           Divorced Handlers-cleaners Not-in-family White
4   3          Private 234721      11th             7 Married-civ-spouse Handlers-cleaners       Husband Black
5   1          Private 338409 Bachelors            13 Married-civ-spouse    Prof-specialty          Wife Black
6   2          Private 284582   Masters            14 Married-civ-spouse   Exec-managerial          Wife White
     sex capitalgain capitalloss hoursperweek native.country class
1   Male           1           0            2  United-States <=50K
2   Male           0           0            0  United-States <=50K
3   Male           0           0            2  United-States <=50K
4   Male           0           0            2  United-States <=50K
5 Female           0           0            2           Cuba <=50K
6 Female           0           0            2  United-States <=50K

Checking class bias now, on column class,

> table(idata$class)

<=50K  >50K 
37155 11687

There is a class bias, a condition observed when the proportion of one event is much smaller than proportion of other event.
Thus, we must sample the observations in approximately equal proportions to get better models.

CREATING TRAINING AND TESTING SAMPLES

#training data
input_ones <- idata[which(idata$class == ">50K"), ]  # all <=50k values
input_zeros <- idata[which(idata$class == "<=50K"), ]  # all <=50k values
set.seed(100)  # for repeatability of samples
input_ones_training_rows <- sample(1:nrow(input_ones), 0.7*nrow(input_ones))  # >50k values for training
input_zeros_training_rows <- sample(1:nrow(input_zeros), 0.7*nrow(input_ones))  # <=50k values training. 
training_ones <- input_ones[input_ones_training_rows, ]  
training_zeros <- input_zeros[input_zeros_training_rows, ]
trainingData <- rbind(training_ones, training_zeros)  # row bind  

# Create Test Data
test_ones <- input_ones[-input_ones_training_rows, ]
test_zeros <- input_zeros[-input_zeros_training_rows, ]
testData <- rbind(test_ones, test_zeros)  # row bind 

One simple optional step is to create Weight of Evidence (WOE) for categorical variables.
This is a very simple, yet useful method - Read more here.

idata_woe <- idata
str(idata_woe)

table(idata_woe$class)

library(car)
library(dplyr)

idata_woe$class <- dplyr::recode(idata_woe$class,"<=50K" = 0,">50K" = 1)
#idata_woe$class <- car::recode(idata_woe$class,"<=50K" = 0,">50K" = 1) - this is not working
idata_woe

i <- 10
q <- quantile(idata_woe$class,
              probs = c(1:(i-1)/i),
              na.rm = TRUE,
              type=3)
cuts <- unique(q)


WOE_Class <- table(findInterval(idata_woe$class,vec = cuts,rightmost.closed = FALSE),idata_woe$class)
WOE_Class <- as.data.frame.matrix(WOE_Class)
WOE_Class$`0`<-rowSums(WOE_Class)
WOE_Class$WOE <- log((WOE_Class$`1`*sum(WOE_Class$`0`))/(WOE_Class$`0`*sum(WOE_Class$`1`)))

> WOE_Class
      0     1      WOE
1 37155     0     -Inf
2 11687 11687 1.430113

Now that we have already calculated WOE, it is only logical that we calculate the Information Value (IV) also,

library(scorecard)

iv(idata_woe, y="class")

> iv(idata_woe, y="class")
          variable info_value
 1:   relationship 1.52771659
 2: marital.status 1.34867700
 3:            age 0.89506214
 4:    capitalgain 0.89314987
 5:     occupation 0.76787908
 6:      education 0.73690518
 7:  education.num 0.73690518
 8:   hoursperweek 0.40824961
 9:         fnlwgt 0.33755308
10:            sex 0.30052143
11:      workclass 0.17056215
12:    capitalloss 0.11592250
13: native.country 0.07509682
14:           race 0.06784214

Now, we have the information value (meaning the strength of relation between “class” and other parameters.
We now start building the logistic regression model,

logitMod <- glm(class ~ relationship + age + capitalgain + occupation + education.num, data=trainingData, family=binomial(link="logit"))
predicted <- plogis(predict(logitMod, testData))  # predicted scores
# or
predicted <- predict(logitMod, testData, type="response")  # predicted scores
summary(logitMod)

Result is as follows,

> summary(logitMod)

Call:
glm(formula = class ~ relationship + age + capitalgain + occupation + 
    education.num, family = binomial(link = "logit"), data = trainingData)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.3583  -0.5408   0.0015   0.6281   3.3139  

Coefficients:
                            Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 -4.43416    0.18263 -24.280  < 2e-16 ***
relationshipNot-in-family   -2.26671    0.05846 -38.776  < 2e-16 ***
relationshipOther-relative  -2.52271    0.19614 -12.862  < 2e-16 ***
relationshipOwn-child       -3.43974    0.14067 -24.452  < 2e-16 ***
relationshipUnmarried       -2.66316    0.09625 -27.670  < 2e-16 ***
relationshipWife             0.31766    0.09130   3.479 0.000503 ***
age                          0.34161    0.02041  16.734  < 2e-16 ***
capitalgain                  0.95230    0.03972  23.974  < 2e-16 ***
occupationAdm-clerical       1.08011    0.13682   7.895 2.91e-15 ***
occupationArmed-Forces       2.89421    1.34843   2.146 0.031844 *  
occupationCraft-repair       1.29746    0.13042   9.949  < 2e-16 ***
occupationExec-managerial    2.08246    0.13025  15.988  < 2e-16 ***
occupationFarming-fishing    0.23595    0.17661   1.336 0.181539    
occupationHandlers-cleaners  0.58579    0.18067   3.242 0.001185 ** 
occupationMachine-op-inspct  0.83323    0.14885   5.598 2.17e-08 ***
occupationOther-service      0.16679    0.16294   1.024 0.306012    
occupationPriv-house-serv   -0.02904    0.74544  -0.039 0.968925    
occupationProf-specialty     1.70451    0.13280  12.836  < 2e-16 ***
occupationProtective-serv    1.72939    0.17889   9.667  < 2e-16 ***
occupationSales              1.56500    0.13284  11.781  < 2e-16 ***
occupationTech-support       1.48003    0.16098   9.194  < 2e-16 ***
occupationTransport-moving   1.21363    0.14650   8.284  < 2e-16 ***
education.num                0.29476    0.01128  26.127  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 22680  on 16359  degrees of freedom
Residual deviance: 13115  on 16337  degrees of freedom
AIC: 13161

Number of Fisher Scoring iterations: 6

Cheers!

Markov Chains - A Mathematical Treatise

Markov Chains

We have a set of stated, $\mathcal{S}$ = {${s_1, s_2, ..., s_m}$}. A Markov process starts in one of these states and moves successively from one state to another. Each move is called a step . If the chain is currently in state $s_i$, then it moves to state $s_j$ at the next step with a probability denoted by $p_{ij}$, and this probability does not depend upon which states the chain was in before the current state.These probabilities $p_{ij}$ are called transition probabilities. The process can remain in the state it is in, and this occurs with probability $p_{ii}$ . Initially a probability is defined by specifying a particular state as the starting state.

Say, we stay in Patna (Capital of the State of Bihar, India), and at my home we usually eat either Roti ($\text{Ro}$), Rice ($\text{Ri}$), Parantha ($\text{Pa}$)) or Bread ($\text{Br}$) for dinner. There are also some days when we eat one particular dinner on two or three or more consecutive days.

Now, we define transition probabilities of using a Transition Matrix which is given as,
$$\text{T}=\begin{array}{c c} & \begin{array}{c c c} \text{Ri} & \text{Ro} & \text{Pa} & \text{Br} \\ \end{array} \\ \begin{array}{c c c} \text{Ri} \\ \text{Ro} \\ \text{Pa} \\ \text{Br}\\ \end{array} & \left[ \begin{array}{c c c} \frac{1}{2} & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \\ 0 & \frac{1}{2} & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{2} & 0 & \frac{1}{4} & \frac{1}{4} \\ \end{array} \right] \end{array}$$
From the above matrix, it can be seen that $P_{(Rice \implies Rice)}=\frac{1}{2}$ and similarly, $P_{(Roti \implies Parantha)}=\frac{1}{4}$.

In the similar lines of the above concept, a Transition Diagram can also be drawn.

Let us dig further into the concept now. Imagine today we eat Rice at home and two days from today we eat Bread. This is written as $P^{2}_{(Rice \implies Bread)}$. We see that if we ate Rice today then the event of us eating Bread two days from today is a disjoint union of the following,

1. We ate Rice tomorrow and Bread day after.
2. We ate something else tomorrow and Bread day after.
3. We ate Bread tomorrow and Bread day after.

Therefore, we can write $P^{2}_{(Rice \implies Bread)}$ as ,
$$P^{2}_{(Rice \implies Bread)} = [P_{(Rice \implies Rice)}]\cdot[P_{(Rice \implies Bread)}]+[P_{(Rice \implies Roti)}]\cdot[P_{(Roti \implies Bread)}]+ \\ [P_{(Rice \implies Parantha)}]\cdot[P_{(Parantha \implies Bread)}]+[P_{(Rice \implies Bread)}]\cdot[P_{(Bread \implies Bread)}]$$
The has pointed towards a more generalized concept of Dot product of two vectors.

Considering there are $r$ states in the Markov chain,
$$P^{2}_{ij}=\sum_{k=1}^{r}P_{ik}P_{kj}$$
This study was for $n=2$ , meaning we are yet to define this in a more generalized manner. Obviously, looking at the present scenario of high-end computing software packages available, this attempt seems to be a waste of efforts, but then knowing the math is always good fun!

Before doing that, we will again go into the basic definition as defined above and try to make an iterative approach into the problem.

We have assumed that our meal starts at state Rice, so, for us,
$$\text{p}_0=\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0\\ \end{bmatrix}$$
$\text{p}_0$ is the state of our system at $0^{th}$ state or the beginning of the system.

Similarly, if we want to find $\text{p}_1$ we will do it as,
$$\text{p}_1 = \text{T} \cdot \text{p}_0$$
In the similar lines,
$$\boxed{\text{p}_5 = \text{T} \cdot \text{p}_4 = \text{T}^{2} \cdot \text{p}_3 = \text{T}^3 \cdot \text{p}_2 = \text{T}^4 \cdot \text{p}_1 = \text{T}^{5} \cdot \text{p}_0}$$
This means, for generalization’s sake,
$$\text{p}_n = \text{T}^{n}\cdot \text{p}_0$$
Using some very basic Linear Algebra methods to compute $\text{T}^n$,

Recall : $\text{T}=\mathcal{u}\cdot D \cdot\mathcal{u}^{-1}$

Here, $D$ is a Diagonal matrix and $\mathcal{u}$ is a matrix whose columns correspond to the Eigen-vectors of $\text{T}$.

I leave the computation to your own practice.

There for,
$$\text{T}^n = \mathcal{u} \cdot D^n . \mathcal{u}^{-1}$$
Since, $D$ is a diagonal matrix, which is obviously of form,
$$D = \begin{bmatrix} a & 0 & 0 & 0 \\ 0 & b & 0 & 0 \\ 0 & 0 & c & 0 \\ 0 & 0 & 0 & d \end{bmatrix}$$
Therefore, $D^n$ will be,
$$D = \begin{bmatrix} a^n & 0 & 0 & 0 \\ 0 & b^n & 0 & 0 \\ 0 & 0 & c^n & 0 \\ 0 & 0 & 0 & d^n \end{bmatrix}$$
Hence, we can easily compute $T^n$.

The Python implementation of this concept is relatively simple once we understand the math behind it.

All the best with that.

Cheers!

Resources:

1. Matrix Diagonalization

2. Eigen Values and Eigen Vectors

3. Markov Chains

4. Markov Chains

Catalan's Conjecture - A learning exercise for the bored mind - Contd.

Before we move to Mihailescu’s proof, let us cover some more mathematical concepts (so that we feel superior to the people around us, just kidding :smile:),

There is an interesting theorem that Mihailescu uses in his proof, called the Stickelberger’s theorem,

---

Stickelberger’s theorem:

This is a result of algebraic number theory, which gives more information about Galois Module structure of class groups of Cyclotomic Fields.
This theorem consists of Stickelberger’s element and Stickelberger’s ideal.
I will now state the complete definition of the theorem and then visit its corners as we move along,

Let $K_m$ denote the $m$-th cyclotomic field${[1]}$ . It is a Galois extension of $\mathbb{Q}$ with Galois Group $G_m$ isomorphic to the multiplicative group of integers modulo $m$ $(\mathbb{Z}/m\mathbb{Z})^{\times}$.

The Stickelberger element (of level $m$ or of $K_m$) is an element in thr group ring $\mathbb{Q}[G_m]$ and the Stickelberger ideal is an ideal in the group ring $\mathbb{Z}[G_m]$. (Note: $\mathbb{Z} \subset \mathbb{Q}$).

The definition of both the Stickelberger element and ideal are, let $\zeta_m$ denote a primitive $m$-th root of unity ${[2]}$. The isomorphism from $(\mathbb{Z/mZ})^{\times}$ to $G_m$ is given by sending $a$ to $\sigma_a$ by the relation $$\sigma_a(\zeta_m)=\zeta^{a}_m.$$
The Stickelberger element of level $m$ is given by,
$$\theta(K_m)=\frac{1}{m}\sum_{a=1 \ (a,m)=1}^ma\cdot \sigma_a^{-1} \in \mathbb{Q}[G_m].$$
The Stickelberger ideal of level $m$ is given by,
$$I(K_m)=\theta(K_m) \mathbb{Z}[G_m] \cap \mathbb{Z}[G_m].$$

---

Inkeri$[3]$, used the concept of Wiefrich pair (explain in the previous blogpost of this series) in the context of Catalan’s equation as follows:

A Wieferich pair is a pair $(p,q)$ of primes such that $$p^{q-1} \equiv 1( \text{mod} \ q^2)$$ and $$q^{p-1} \equiv 1(\text{mod} \ p^2).$$

He showed that if the Catalan’s equation holds, then either $(p,q)$ is a Wieferich pair, or $q$ divides $h_p$, the class number$[4]$ of cyclotomic field $\mathbb{Q}(\zeta_p)$, or $p$ divides $h_q$, the class number of cyclotomic field $\mathbb{Q}(\zeta_q)$, there were other developments in this direction too.

Bugeaud and Hanrot $[5]$proved a class number criterion concerning Catalan’s equation, which implies that the Catalan’s Equation ( $x^p - y^q =1$) has no solutions in non-zero integers $x$ and $y$ if $p$ and $q$ are primes such that one of them is smaller than 43. This was a huge achievement , I recommend you to have a look at the paper at [5].

Mihailescu proved that the Catalan equation has no solutions if $p$ and $q$ are odd and $q$ does not divide $p - 1$. By this result the Catalan conjecture became a theorem. And later Mihailescu succeeded in finding a more elegant proof of Catalan’s conjecture in the case where $q$ does divide $p-1$. Thus, Catalan’s conjecture is a theorem with an algebraic proof in which no computer calculations
are needed.

---

---

In this section e wll discuss some breakthrough results by Mihailescu. The most important one is that $q^2$ divides $x$.

The following lemma will be used for that, an element $a$ of a ring $R$ is called nilpotent if $\exists$ and integer $n$ such that $a^n=0$.

Lemma 1: The ring $\mathbb{O}_k/(q)$ does not contain nilpotent elements, if $\alpha$, $\beta$ $\in$ $\mathbb{O}_k = \mathbb{Z}[\zeta]$ satisfy the congruence $\alpha^p \equiv \beta^q (\text{mod} \ q^2)$

Theorem 1: For $\theta \in {I_s}^{-}$, the element $(x-\zeta)^\theta$ is a $q^{th}$ in $K=\mathbb{Q}(\zeta)$. We also have that $q^2$ divides $x$ and $p^2$ divides $y$.

The proofs of the above Lemma and Theorem are beyond the scope of this blog; regardless to say that the pre-requisites are already covered in detail. For the interested, you can refer Catalan’s Conjecture - A cyclotomic field.

Cheers!

Catalan's Conjecture - A learning exercise for the bored mind.

I was going through a video by Numberphile where they were talking about the Catalan’s conjecture.
This is another conjecture which is simple to state; however extremely difficult to prove; just like Collatz conjecture, about which I already have a blog in place.
At the very outset of the blog; let me bring it to your knowledge that this blog is just for learning and understanding and contains “very little”-to-“null” original work on the subject. However, this blog will be exhaustive and will contain a lot of references that will help a newbie (like myself) get into the depths of the conjecture and fully understand concepts used in its proof.
Let us now define the conjecture;
This statement was conjectured by Eugene Catalan (1814–1894) and was sent to the editor of Journal fur die Reine und Angewandte Mathematik$^{[1]-\text{Introduction}}$.
The conjecture is as follows,
$2^3$ and $3^2$are two powers of natural numbers whose values are consecutive (i.e., 8 and 9); the conjecture is for a mathematical statement such as,
$$\boxed{x^a-y^b=1} \tag*{(1)}$$
The only solution of $(1)$ is for $a$, $b$ $>$ $1$ and $x$, $y$ $>$ $0$ is $x=3$, $a=2$, $y=2$, $b=3$.
Catalan’s conjecture was proven true by Preda Mihăilescu in 2002; the proof involves the theory of cyclotomic fields and Galois Modules$^{[2]-\text{History}}$.
So, as you see now, the breadth of the subject blew up! On second were talking about squares and cubes and now we are talking of Cyclotomic fields and Galois Modules!
Nevertheless, I will try and cover each topic in brief and quench our mathematical thirst!
Let us consider $(1)$ (with $a=p$ and $b=q$; for the sake of eliminating any confusion between $a$ and an english “a”), unless otherwise stated, $x$ and $y$ can be negative integers as well. Now, we re-write $(1)$ as follows,
$$(x-1)\frac{x^p-1}{x-1}=y^q \tag*{(2)}$$
The GCD of the two factors on the left hand side of the equation (after considering $x^p = ((x-1)+1)^p$) is either $p$ or $1$ (How did this happen?).

---

Some concepts before we move ahead:

The Wieferich pairs
$^{[3]-\text{Wieferich pair}}$In mathematics,
a Wieferich pair is a pair of prime number $p$ and $q$ that satisfy,
$$p^{q-1} \equiv 1( \mod p^2)$$
Let us write $(2)$ as, $$(x-1)\frac{x^p-1}{x-1}=pb^q$$ for, $p \nmid b$
This suggests a traditional approach of factorizing the left hand side
in $\mathbb{Z}[\zeta]$, the ring of integers in the $p$th cyclotomic
field $\mathbb{Q}(\zeta)(\zeta=e^{2\pi i/p})$

---

Ring$^{[4]-\text{Ring-Wolfram}}$:
A ring in the mathematical sense is a set $S$ together with two binary
operators $+$ and $.$ (Addition and multiplication), satisfying the
following conditions:
1. Additive associativity: For all $a,b,c \in S$, $(a+b)+c=(a+b)+c$,
2. Additive commutativity: For all $a,b\in S$, $a+b=b+a$,
3. Additive identity: There exists an element $0$ in $S$ such that for all a in $S$, $0+a=a+0=a$.
4. Additive inverse: For every a in $S$ there exists $-a \in S$ such that $a+(-a)=(-a)+a=0$,
5. Left and right distributivity: For all $a,b,c \in S$, $a.(b+c)=(a.b)+(a.c)$ and $(b+c).a=(b.a)+(c.a)$,
6. Multiplicative associativity: For all $a,b,c \in S$, $(a.b).c=a.(b.c)$ ($a$ ring satisfying this property is sometimes
explicitly termed an associative ring). Conditions 1-5 are always
required. Though non-associative rings exist, virtually all texts also
require condition 6.
7. Multiplicative commutativity: For all $a,b \in S$, $a.b=b.a$ ($a$ ring satisfying this property is termed a commutative ring),
8. Multiplicative identity: There exists an element $1 \in S$ such that for all $a \neq 0 \in S$, $1.a=a.1=a$ (a ring satisfying this
property is termed a unit ring, or sometimes a “ring with identity”),
9. Multiplicative inverse: For each $a \neq 0 \in S$, there exists an element $a^{-1} \in S$ such that $\forall a \neq 0 \in S, a.a^{-1}=a^{-1}.a=1$, where $1$ is the identity element.

---

A brief history of the past developments on this conjecture is a must to be read and understood; the significance comes due to the period of 150 years for which it remained an open problem$^{[5]-\text{Catalan's Conjecture: Are 8 and 9 the Only Consecutive Powers?}}$,

• Only after six years after Catalan formally defined the conjecture, a result was proposed by French mathematician, Victor Lebesgue. He stated that, for the equation, $x^p-y^2=1$; where $p$ is a prime; has no solutions for positive values of $x$ and $y$. A proof of the same will be discussed in brief in the later part of the blog.
• After Lebesgue’s work, all development solely consisted of small exponents, and then Naggel showed in 1921 that the difference between a third power and an other perfect power never is equal to 1.
• In 1932, Selberg proved that, $x^4-y^n=1$ has no solution in positive integers when $n>1$. A stronger result to this was proved by Ko Cho in 1965, that stated that the equation $x^2-y^q=1$ has no solutions for positive integers when $q \geq 5$.
• Cassels made some observations for $x^p-y^q=1$ where $p$ and $q$ are odd-primes. He proved that is this equality holds for positive integers $x$ and $y$, then $p$ divides $y$ and $q$ divides $x$. For the case $p=2$, this had already been shown by Naggel
• Inkeri defined the concept of a Wieferich pair [the definition and explanation of the same is given above] in the concept of Catalan equation as follows:
  If the Catalan’s equation $(1)$ holds, then either $(p,q)$ is a Wieferich Pair, or $q$ divides $h_p$, the class number of the cyclotomic field $\mathbb{Q}(\zeta_p)$, or $p$ divides $h_q$, the class number of the cyclotomic field $\mathbb{Q}(\zeta_q)$

---

Cyclotomic Field$^{[5]-\text{Cyclotomic field}}$:
In number theory, a cyclotomic field is a number field obtained by adjoining a complex primitive root of unity to $\mathbb{Q}$, the field of rational numbers. The $n$-th cyclotomic field $\mathbb{Q}(\zeta_n)$ (where $n > 2$) is obtained by adjoining a primitive $n$-th root of
Primitive root of unity$^{[6]-\text{Root of unity}}$
In mathematics, a root of unity, occasionally called a de Moivre number, is any complex number that gives $1$ when raised to some positive integer power $n$.

---

• Some time later Mihailescu proved that the Catalan equation has no solutions if $p$ and $q$ are odd and $q$ does not divide $p-1$. By this result the Catalan conjecture became a theorem

---

---

More on Cyclotomic Fields:
Let $p$ be an odd-prime number. Let $\Phi_p$ be the $p$-th cyclotomic polynomial in $\mathbb{Q}[X]$ i.e., $\Phi_p = \frac{X^p-1}{X-1}$. Consider the field extension $\mathbb{Q}[X]/(\Phi_p) \cong \mathbb{Q}(\zeta)$ of $\mathbb{Q}$, where $\zeta$ denotes a primitive $p$-th root of unity. This is a field extension of degree $p-1$ and it is reducible in $\mathbb{Q}[X]$. We denote $\mathbb{Q}(\zeta)$ by $k$ from now on.
This field extension is Galois with Galois Group,
$$G=\text{Gal}(\mathbb{Q}(\zeta)/\mathbb{Q}) \cong (\mathbb{Z}/p\mathbb{Z})^{*},$$
Since the map,
$$(\mathbb{Z}/p\mathbb{Z})^{*} {{\sim} \atop {\to}}\text{Gal}(\mathbb{Q}(\zeta)/\mathbb{Q})$$
$$a (\mod p) \mapsto (\sigma_a : \zeta \mapsto \zeta^a)$$
is an isomorphism
The automorphism $\sigma_{p-1}$ acts in all embeddings as complex conjugation. Therefore, we call $\sigma_{p-1}$ complex conjugation.
The fixed field of complex conjugation is $\mathbb{Q}(\zeta+\zeta^{-1})$, which is called the maximal real subfield of $\mathbb{Q}(\zeta)$. We denote $\mathbb{Q}(\zeta + \zeta^{-1})$ by $K^{+}$. The field extension $\mathbb{Q}(\zeta+\zeta^{-1})$ of $\mathbb{Q}$ has degree $\frac{p-1}{2}$ and it is Galois with Galois theory.
$$G^{+}=\text{Gal}(\mathbb{Q}(\zeta+\zeta^{-1})/\mathbb{Q}) \cong (\mathbb{Z}/p\mathbb{Z})^{*}/(\pm 1)$$

Another important concept that is

Mihailescu’s proof

To be contd…in the next blogpost!

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!