Linear Algebra behind the lm() function in R

This post comes out of the blue, nearly 2 years since my last one. I realize I’ve been lazy, so here’s hoping I move from an inertia of rest to that of motion, implying, regular and (hopefully) relevant posts. I also chanced upon some wisdom while scrolling through my Twitter feed:

This blog post in particular was meant to be a reminder to myself and other R users that the much used lm() function in R (for fitting linear models) can be replaced with some handy matrix operations to obtain regression coefficients, their standard errors and other goodness-of-fit stats printed out when summary() is called on an lm object.

Linear regression can be formulated mathematically as follows:
\mathbf{y} = \mathbf{X} \mathbf{\beta} + \mathbf{\epsilon} ,
\mathbf{\epsilon} \sim N(0, \sigma^2 \mathbf{I})

\mathbf{y} is the \mathbf{n}\times \mathbf{1} outcome variable and \mathbf{X} is the \mathbf{n}\times \mathbf{(\mathbf{k}+1)} data matrix of independent predictor variables (including a vector of ones corresponding to the intercept). The ordinary least squares (OLS) estimate for the vector of coefficients \mathbf{\beta} is:

\hat{\mathbf{\beta}} = (\mathbf{X}^{\prime} \mathbf{X})^{-1} \mathbf{X}^{\prime} \mathbf{y}

The covariance matrix can be obtained with some handy matrix operations:
\textrm{Var}(\hat{\mathbf{\beta}}) = (\mathbf{X}^{\prime} \mathbf{X})^{-1} \mathbf{X}^{\prime} \;\sigma^2 \mathbf{I} \; \mathbf{X} (\mathbf{X}^{\prime} \mathbf{X})^{-1} = \sigma^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}
given that \textrm{Var}(AX) = A \times \textrm{Var}X \times A^{\prime}; \textrm{Var}(\mathbf{y}) = \mathbf{\sigma^2}

The standard errors of the coefficients are basically \textrm{Diag}(\sqrt{\textrm{Var}(\hat{\mathbf{\beta}})}) = \textrm{Diag}(\sqrt{\sigma^2 (\mathbf{X}^{\prime} \mathbf{X})^{-1}}) and with these, one can compute the t-statistics and their corresponding p-values.

Lastly, the F-statistic and its corresponding p-value can be calculated after computing the two residual sum of squares (RSS) statistics:

  • \mathbf{RSS} – for the full model with all predictors
  • \mathbf{RSS_0} – for the partial model (\mathbf{y} = \mathbf{\mu} + \mathbf{\nu}; \mathbf{\mu} = \mathop{\mathbb{E}}[\mathbf{y}]; \mathbf{\nu} \sim N(0, \sigma_0^2 \mathbf{I}) ) with the outcome observed mean as estimated outcome

\mathbf{F} = \frac{(\mathbf{RSS_0}-\mathbf{RSS})/\mathbf{k}}{\mathbf{RSS}/(\mathbf{n}-\mathbf{k}-1)}

I wrote some R code to construct the output from summarizing lm objects, using all the math spewed thus far. The data used for this exercise is available in R, and comprises of standardized fertility measures and socio-economic indicators for each of 47 French-speaking provinces of Switzerland from 1888. Try it out and see for yourself the linear algebra behind linear regression.

### Linear Regression Using lm() ----------------------------------------
dat <- swiss
linear_model <- lm(Fertility ~ ., data = dat)
# Call:
# lm(formula = Fertility ~ ., data = dat)
# Residuals:
# Min 1Q Median 3Q Max
# -15.2743 -5.2617 0.5032 4.1198 15.3213
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
# Agriculture -0.17211 0.07030 -2.448 0.01873 *
# Examination -0.25801 0.25388 -1.016 0.31546
# Education -0.87094 0.18303 -4.758 2.43e-05 ***
# Catholic 0.10412 0.03526 2.953 0.00519 **
# Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Residual standard error: 7.165 on 41 degrees of freedom
# Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
# F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10
### Using Linear Algebra ------------------------------------------------
y <- matrix(dat$Fertility, nrow = nrow(dat))
X <- cbind(1, as.matrix(x = dat[,-1]))
colnames(X)[1] <- "(Intercept)"
# N x k matrix
N <- nrow(X)
k <- ncol(X) - 1 # number of predictor variables (ergo, excluding Intercept column)
# Estimated Regression Coefficients
beta_hat <- solve(t(X)%*%X)%*%(t(X)%*%y)
# Variance of outcome variable = Variance of residuals
sigma_sq <- residual_variance <- (N-k-1)^-1 * sum((y - X %*% beta_hat)^2)
residual_std_error <- sqrt(residual_variance)
# Variance and Std. Error of estimated coefficients of the linear model
var_betaHat <- sigma_sq * solve(t(X) %*% X)
coeff_std_errors <- sqrt(diag(var_betaHat))
# t values of estimates are ratio of estimated coefficients to std. errors
t_values <- beta_hat / coeff_std_errors
# p-values of t-statistics of estimated coefficeints
p_values_tstat <- 2 * pt(abs(t_values), N-k, lower.tail = FALSE)
# assigning R's significance codes to obtained p-values
signif_codes_match <- function(x){
ifelse(x <= 0.001,"***",
ifelse(x <= 0.01,"**",
ifelse(x < 0.05,"*",
ifelse(x < 0.1,"."," "))))
signif_codes <- sapply(p_values_tstat, signif_codes_match)
# R-squared and Adjusted R-squared (refer any econometrics / statistics textbook)
R_sq <- 1 - (N-k-1)*residual_variance / (N*mean((y - mean(y))^2))
R_sq_adj <- 1 - residual_variance / ((N/(N-1))*mean((y - mean(y))^2))
# Residual sum of squares (RSS) for the full model
RSS <- (N-k-1)*residual_variance
# RSS for the partial model with only intercept (equal to mean), ergo, TSS
RSS0 <- TSS <- sum((y - mean(y))^2)
# F statistic based on RSS for full and partial models
# k = degress of freedom of partial model
# N - k - 1 = degress of freedom of full model
F_stat <- ((RSS0 - RSS)/k) / (RSS/(N-k-1))
# p-values of the F statistic
p_value_F_stat <- pf(F_stat, df1 = k, df2 = N-k-1, lower.tail = FALSE)
# stitch the main results toghether
lm_results <-, coeff_std_errors,
t_values, p_values_tstat, signif_codes))
colnames(lm_results) <- c("Estimate","Std. Error","t value","Pr(>|t|)","")
### Print out results of all relevant calcualtions -----------------------
cat("Residual standard error: ",
round(residual_std_error, digits = 3),
" on ",N-k-1," degrees of freedom",
"\nMultiple R-squared: ",R_sq," Adjusted R-squared: ",R_sq_adj,
"\nF-statistic: ",F_stat, " on ",k-1," and ",N-k-1,
" DF, p-value: ", p_value_F_stat,"\n")
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 66.9151816789654 10.7060375853301 6.25022854119771 1.73336561301153e-07 ***
# Agriculture -0.172113970941457 0.0703039231786469 -2.44814177018405 0.0186186100433133 *
# Examination -0.258008239834722 0.253878200892098 -1.01626779663678 0.315320687313066
# Education -0.870940062939429 0.183028601571259 -4.75849159892283 2.3228265226988e-05 ***
# Catholic 0.104115330743766 0.035257852536169 2.95296858017545 0.00513556154915653 **
# Infant.Mortality 1.07704814069103 0.381719650858061 2.82156849475775 0.00726899472564356 **
# Residual standard error: 7.165 on 41 degrees of freedom
# Multiple R-squared: 0.706735 Adjusted R-squared: 0.670971
# F-statistic: 19.76106 on 4 and 41 DF, p-value: 5.593799e-10

Hope this was useful and worth your time!


MITx: 6.008.1x Computational Probability and Inference

I got really interested in Computational Probability and Inference (6.008.1x) for the following reasons:

  1. I love probability and have solved countless problems on probability ever since I learned math
  2. …and yet I’ve never coded up probabilistic models!
  3. The assignments and project work for this course are to be implemented in Python!

You don’t need to have prior experience in either probability or inference, but you should be comfortable with basic Python programming and calculus.

– Basic discrete probability theory
– Graphical models as a data structure for representing probability distributions
– Algorithms for prediction and inference
– How to model real-world problems in terms of probabilistic inference

The course started on September 12, is 12-weeks long and is structured in the following manner:

Week 1 (9/12 – 9/16): Introduction to probability and computation
A first look at basic discrete probability, how to interpret it, what probability spaces and random variables are, and how to code these up and do basic simulations and visualizations.

Week 2 (9/19 – 9/23): Incorporating observations
Incorporating observations using jointly distributed random variables and using events. Three classic probability puzzles are presented to help elucidate how to interpret probability: Simpson’s paradox, Monty Hall, boy or girl paradox.

Week 3 (9/26 – 9/30): Introduction to inference, structure in distributions, and information measures
The product rule and inference with Bayes’ theorem. Independence: A structure in distributions. Measures of randomness: entropy and information divergence. Mutual information.

Week 4 (10/3 – 10/7): Expectations, and driving to infinity in modeling uncertainty
Expected values of random variables. Classic puzzle: the two envelope problem. Probability spaces and random variables that take on a countably infinite number of values and inference with these random variables.

Week 5 (10/10 – 10/14): Efficient representations of probability distributions on a computer
Introduction to undirected graphical models as a data structure for representing probability distributions and the benefits/drawbacks of these graphical models. Incorporating observations with graphical models.

Week 6 (10/17 – 10/21): Inference with graphical models, part I
Computing marginal distributions with graphical models in undirected graphical models including hidden Markov models..

Week 7 (10/24 – 10/28): Inference with graphical models, part II
Computing most probable configurations with graphical models including hidden Markov models.

Week 8 (10/31 – 11/4): Introduction to learning probability distributions
Learning an underlying unknown probability distribution from observations using maximum likelihood. Three examples: estimating the bias of a coin, the German tank problem, and email spam detection.

Week 9 (11/7 – 11/11): Parameter estimation in graphical models
Given the graph structure of an undirected graphical model, we examine how to estimate all the tables associated with the graphical model.

Week 10 (11/14 – 11/18): Model selection with information theory
Learning both the graph structure and the tables of an undirected graphical model with the help of information theory. Mutual information of random variables.

Week 11 (11/21 – 11/25): Final project
Final project assigned

Week 12 (11/28 – 12/2): Final project


I’m SO taking this course. Hope this interests you as well!

Implementing Undirected Graphs in Python

There are 2 popular ways of representing an undirected graph.

Adjacency List
Each list describes the set of neighbors of a vertex in the graph.


Adjacency Matrix
The elements of the matrix indicate whether pairs of vertices are adjacent or not in the graph.


Here’s an implementation of the above in Python:

class Vertex:
def __init__(self, vertex): = vertex
self.neighbors = []
def add_neighbor(self, neighbor):
if isinstance(neighbor, Vertex):
if not in self.neighbors:
self.neighbors = sorted(self.neighbors)
neighbor.neighbors = sorted(neighbor.neighbors)
return False
def add_neighbors(self, neighbors):
for neighbor in neighbors:
if isinstance(neighbor, Vertex):
if not in self.neighbors:
self.neighbors = sorted(self.neighbors)
neighbor.neighbors = sorted(neighbor.neighbors)
return False
def __repr__(self):
return str(self.neighbors)
class Graph:
def __init__(self):
self.vertices = {}
def add_vertex(self, vertex):
if isinstance(vertex, Vertex):
self.vertices[] = vertex.neighbors
def add_vertices(self, vertices):
for vertex in vertices:
if isinstance(vertex, Vertex):
self.vertices[] = vertex.neighbors
def add_edge(self, vertex_from, vertex_to):
if isinstance(vertex_from, Vertex) and isinstance(vertex_to, Vertex):
if isinstance(vertex_from, Vertex) and isinstance(vertex_to, Vertex):
self.vertices[] = vertex_from.neighbors
self.vertices[] = vertex_to.neighbors
def add_edges(self, edges):
for edge in edges:
def adjacencyList(self):
if len(self.vertices) >= 1:
return [str(key) + ":" + str(self.vertices[key]) for key in self.vertices.keys()]
return dict()
def adjacencyMatrix(self):
if len(self.vertices) >= 1:
self.vertex_names = sorted(g.vertices.keys())
self.vertex_indices = dict(zip(self.vertex_names, range(len(self.vertex_names))))
import numpy as np
self.adjacency_matrix = np.zeros(shape=(len(self.vertices),len(self.vertices)))
for i in range(len(self.vertex_names)):
for j in range(i, len(self.vertices)):
for el in g.vertices[self.vertex_names[i]]:
j = g.vertex_indices[el]
self.adjacency_matrix[i,j] = 1
return self.adjacency_matrix
return dict()
def graph(g):
""" Function to print a graph as adjacency list and adjacency matrix. """
return str(g.adjacencyList()) + '\n' + '\n' + str(g.adjacencyMatrix())
a = Vertex('A')
b = Vertex('B')
c = Vertex('C')
d = Vertex('D')
e = Vertex('E')
g = Graph()


["A:['B', 'C', 'E']", "C:['A', 'B', 'D', 'E']", "B:['A', 'C', 'D']", "E:['A', 'C']", "D:['B', 'C']"]
[[ 0. 1. 1. 0. 1.]
[ 1. 0. 1. 1. 0.]
[ 1. 1. 0. 1. 1.]
[ 0. 1. 1. 0. 0.]
[ 1. 0. 1. 0. 0.]]

Computing Work Done (Total Pivot Comparisons) by Quick Sort

A key aspect of the Quick Sort algorithm is how the pivot element is chosen. In my earlier post on the Python code for Quick Sort, my implementation takes the first element of the unsorted array as the pivot element.

However with some mathematical analysis it can be seen that such an implementation is O(n2) in complexity while if a pivot is randomly chosen, the Quick Sort algorithm is O(nlog2n).

To witness this in action, one can measure the work done by the algorithm comparing two cases, one with a randomized pivot choice – and one with a fixed pivot choice, say the first element of the array (or the last element of the array).


A decent proxy for the amount of work done by the algorithm would be the number of pivot comparisons. These comparisons needn’t be computed one-by-one, rather when there is a recursive call on a subarray of length m, you should simply add m−1 to your running total of comparisons.

3 Cases

To put things in perspective, let’s look at 3 cases. (This is basically straight out of a homework assignment from Tim Roughgarden’s course on the Design and Analysis of Algorithms).
Case I with the pivot being the first element.
Case II with the pivot being the last element.
Case III using the “median-of-three” pivot rule. The primary motivation behind this rule is to do a little bit of extra work to get much better performance on input arrays that are nearly sorted or reverse sorted.

Median-of-Three Pivot Rule

Consider the first, middle, and final elements of the given array. (If the array has odd length it should be clear what the “middle” element is; for an array with even length 2k, use the kth element as the “middle” element. So for the array 4 5 6 7, the “middle” element is the second one —- 5 and not 6! Identify which of these three elements is the median (i.e., the one whose value is in between the other two), and use this as your pivot.

Python Code

This file contains all of the integers between 1 and 10,000 (inclusive, with no repeats) in unsorted order. The integer in the ith row of the file gives you the ith entry of an input array. I downloaded this file and named it QuickSort_List.txt

You can run the code below and see for yourself that the number of comparisons for Case III are 138,382 compared to 162,085 and 164,123 for Case I and Case II respectively. You can play around with the code in an IPython / Jupyter notebook here.

# Case I
# First element of the unsorted array is chosen as pivot element for sorting using Quick Sort
def countComparisonsWithFirst(x):
""" Counts number of comparisons while using Quick Sort with first element of unsorted array as pivot """
global count_pivot_first
if len(x) == 1 or len(x) == 0:
return x
count_pivot_first += len(x)-1
i = 0
for j in range(len(x)-1):
if x[j+1] < x[0]:
x[j+1],x[i+1] = x[i+1], x[j+1]
i += 1
x[0],x[i] = x[i],x[0]
first_part = countComparisonsWithFirst(x[:i])
second_part = countComparisonsWithFirst(x[i+1:])
return first_part + second_part
# Case II
# Last element of the unsorted array is chosen as pivot element for sorting using Quick Sort
def countComparisonsWithLast(x):
""" Counts number of comparisons while using Quick Sort with last element of unsorted array as pivot """
global count_pivot_last
if len(x) == 1 or len(x) == 0:
return x
count_pivot_last += len(x)-1
x[0],x[-1] = x[-1],x[0]
i = 0
for j in range(len(x)-1):
if x[j+1] < x[0]:
x[j+1],x[i+1] = x[i+1], x[j+1]
i += 1
x[0],x[i] = x[i],x[0]
first_part = countComparisonsWithLast(x[:i])
second_part = countComparisonsWithLast(x[i+1:])
return first_part + second_part
# Case III
# Median-of-three method used to choose pivot element for sorting using Quick Sort
def middle_index(x):
""" Returns the index of the middle element of an array """
if len(x) % 2 == 0:
middle_index = len(x)/2 - 1
middle_index = len(x)/2
return middle_index
def median_index(x,i,j,k):
""" Returns the median index of three when passed an array and indices of any 3 elements of that array """
if (x[i]-x[j])*(x[i]-x[k]) < 0:
return i
elif (x[j]-x[i])*(x[j]-x[k]) < 0:
return j
return k
def countComparisonsMedianOfThree(x):
""" Counts number of comparisons while using Quick Sort with median-of-three element is chosen as pivot """
global count_pivot_median
if len(x) == 1 or len(x) == 0:
return x
count_pivot_median += len(x)-1
k = median_index(x, 0, middle_index(x), -1)
if k != 0: x[0], x[k] = x[k], x[0]
i = 0
for j in range(len(x)-1):
if x[j+1] < x[0]:
x[j+1],x[i+1] = x[i+1], x[j+1]
i += 1
x[0],x[i] = x[i],x[0]
first_part = countComparisonsMedianOfThree(x[:i])
second_part = countComparisonsMedianOfThree(x[i+1:])
return first_part + second_part
# initializing counts
count_pivot_first = 0; count_pivot_last = 0; count_pivot_median = 0
# Cast I
# Read the contents of the file into a Python list
NUMLIST_FILENAME = "QuickSort_List.txt"
inFile = open(NUMLIST_FILENAME, 'r')
with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
# call functions to count comparisons
# Read the contents of the file into a Python list
NUMLIST_FILENAME = "QuickSort_List.txt"
inFile = open(NUMLIST_FILENAME, 'r')
with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
# call functions to count comparisons
# Read the contents of the file into a Python list
NUMLIST_FILENAME = "QuickSort_List.txt"
inFile = open(NUMLIST_FILENAME, 'r')
with inFile as f: numList = [int(integers.strip()) for integers in f.readlines()]
# call functions to count comparisons
print count_pivot_first, count_pivot_last, count_pivot_median

Generating Permutation Matrices in Octave / Matlab

I have been doing Gilbert Strang’s linear algebra assignments, some of which require you to write short scripts in MatLab, though I use GNU Octave (which is kind of like a free MatLab). I was trying out this problem:

permutationMatricesTo solve this quickly, it would have been nice to have a function that would give a list of permutation matrices for every n-sized square matrix, but there was none in Octave, so I wrote a function permMatrices which creates a list of permutation matrices for a square matrix of size n.

% function to generate permutation matrices given the size of the desired permutation matrices
function x = permMatrices(n)
x = zeros(n,n,factorial(n));
permutations = perms(1:n);
for i = 1:size(x,3)
x(:,:,i) = eye(n)(permutations(i,:),:);
view raw permMatrices.m hosted with ❤ by GitHub

For example:


The MatLab / Octave code to solve this problem is shown below:

% Solution for part (a)
p = permMatrices(3);
n = size(p,3); % number of permutation matrices
v = zeros(n,1); % vector of zeros with dimension equalling number of permutation matrices
% check for permutation matrices other than identity matrix with 3rd power equalling identity matrix
for i = 1:n
if p(:,:,i)^3 == eye(3)
v(i,1) = 1;
v(1,1) = 0; % exclude identity matrix
ans1 = p(:,:,v == 1)
% Solution for part (b)
P = permMatrices(4);
m = size(P,3); % number of permutation matrices
t = zeros(m,1); % vector of zeros with dimension equalling number of permutation matrices
% check for permutation matrices with 4th power equalling identity matrix
for i = 1:m
if P(:,:,i)^4 == eye(4)
t(i,1) = 1;
% print the permutation matrices
ans2 = P(:,:,t == 0)
view raw Section2_7_13.m hosted with ❤ by GitHub


Output for 13(a)
Output for 13(b)


Sherlock and the Beast – HackerRank

I found myself stuck on this problem recently. I must confess, I lost a couple of hours trying to get to figure the logic for this one. Here’s the problem:


I’ve written 2 functions to solve this problem. The first one I used for smaller N, say N < 30 and the second one for N > 30. The second function is elegant, and it relies on the mathematical property that if a number N is not divisible by 3, it could either leave a remainder 1 or 2.

If it leaves a remainder 2, then subtracting 5 once would make the number divisible by 3. If it leaves a remainder 1, then subtracting 5 twice would make the number divisible by 3.

We subtract 5 from N iteratively and attempt to divide N into 2 parts, one divisible by 3 and the other divisible by 5. We want the part that is divisible by 3 to be the larger part, so that the associated Decent Number is the largest possible. This explanation might seem obtuse, but if you get pen down on paper, you’ll understand what I mean.



Karatsuba Multiplication Algorithm – Python Code

Motivation for this blog post

I’ve enrolled in Stanford Professor Tim Roughgarden’s Coursera MOOC on the design and analysis of algorithms, and while he covers the theory and intuition behind the algorithms in a surprising amount of detail, we’re left to implement them in a programming language of our choice.

And I’m ging to post Python code for all the algorithms covered during the course!

The Karatsuba Multiplication Algorithm

Karatsuba’s algorithm reduces the multiplication of two n-digit numbers to at most  n^{\log_23}\approx n^{1.585} single-digit multiplications in general (and exactly n^{\log_23} when n is a power of 2). Although the familiar grade school algorithm for multiplying numbers is how we work through multiplication in our day-to-day lives, it’s slower (\Theta(n^2)\,\!) in comparison, but only on a computer, of course!

Here’s how the grade school algorithm looks:
(The following slides have been taken from Tim Roughgarden’s notes. They serve as a good illustration. I hope he doesn’t mind my sharing them.)


…and this is how Karatsuba Multiplication works on the same problem:



A More General Treatment

Let x and y be represented as n-digit strings in some base B. For any positive integer m less than n, one can write the two given numbers as

x = x_1B^m + x_0
y = y_1B^m + y_0,

where x_0 and y_0 are less than B^m. The product is then

xy = (x_1B^m + x_0)(y_1B^m + y_0)
xy = z_2B^{2m} + z_1B^m + z_0


z_2 = x_1y_1
z_1 = x_1y_0 + x_0y_1
z_0 = x_0y_0

These formulae require four multiplications, and were known to Charles Babbage. Karatsuba observed that xy can be computed in only three multiplications, at the cost of a few extra additions. With z_0 and z_2 as before we can calculate

z_1 = (x_1 + x_0)(y_1 + y_0) - z_2 - z_0

which holds since

z_1 = x_1y_0 + x_0y_1
z_1 = (x_1 + x_0)(y_1 + y_0) - x_1y_1 - x_0y_0

A more efficient implementation of Karatsuba multiplication can be set as xy = (b^2 + b)x_1y_1 - b(x_1 - x_0)(y_1 - y_0) + (b + 1)x_0y_0, where b = B^m.


To compute the product of 12345 and 6789, choose B = 10 and m = 3. Then we decompose the input operands using the resulting base (Bm = 1000), as:

12345 = 12 · 1000 + 345
6789 = 6 · 1000 + 789

Only three multiplications, which operate on smaller integers, are used to compute three partial results:

z2 = 12 × 6 = 72
z0 = 345 × 789 = 272205
z1 = (12 + 345) × (6 + 789) − z2z0 = 357 × 795 − 72 − 272205 = 283815 − 72 − 272205 = 11538

We get the result by just adding these three partial results, shifted accordingly (and then taking carries into account by decomposing these three inputs in base 1000 like for the input operands):

result = z2 · B2m + z1 · Bm + z0, i.e.
result = 72 · 10002 + 11538 · 1000 + 272205 = 83810205.

Pseudocode and Python code

procedure karatsuba(num1, num2)
if (num1 < 10) or (num2 < 10)
return num1*num2
/* calculates the size of the numbers */
m = max(size_base10(num1), size_base10(num2))
m2 = m/2
/* split the digit sequences about the middle */
high1, low1 = split_at(num1, m2)
high2, low2 = split_at(num2, m2)
/* 3 calls made to numbers approximately half the size */
z0 = karatsuba(low1,low2)
z1 = karatsuba((low1+high1),(low2+high2))
z2 = karatsuba(high1,high2)
return (z2*10^(2*m2))+((z1-z2-z0)*10^(m2))+(z0)

def karatsuba(x,y):
"""Function to multiply 2 numbers in a more efficient manner than the grade school algorithm"""
if len(str(x)) == 1 or len(str(y)) == 1:
return x*y
n = max(len(str(x)),len(str(y)))
nby2 = n / 2
a = x / 10**(nby2)
b = x % 10**(nby2)
c = y / 10**(nby2)
d = y % 10**(nby2)
ac = karatsuba(a,c)
bd = karatsuba(b,d)
ad_plus_bc = karatsuba(a+b,c+d) - ac - bd
# this little trick, writing n as 2*nby2 takes care of both even and odd n
prod = ac * 10**(2*nby2) + (ad_plus_bc * 10**nby2) + bd
return prod
view raw hosted with ❤ by GitHub

Magic 5-gon Ring — Project Euler (Problem 68)

Yet another exciting math problem that requires an algorithmic approach to arrive at a quick solution! There is a pen-paper approach to it too, but this post assumes we’re more interested in discussing the programming angle.

First, the problem:

Working clockwise, and starting from the group of three with the numerically lowest external node (4,3,2 in this example), each solution can be described uniquely. For example, the above solution can be described by the set: 4,3,2; 6,2,1; 5,1,3.

It is possible to complete the ring with four different totals: 9, 10, 11, and 12. There are eight solutions in total.

Total Solution Set:
9 4,2,3; 5,3,1; 6,1,2
9 4,3,2; 6,2,1; 5,1,3
10 2,3,5; 4,5,1; 6,1,3
10 2,5,3; 6,3,1; 4,1,5
11 1,4,6; 3,6,2; 5,2,4
11 1,6,4; 5,4,2; 3,2,6
12 1,5,6; 2,6,4; 3,4,5
12 1,6,5; 3,5,4; 2,4,6
By concatenating each group it is possible to form 9-digit strings; the maximum string for a 3-gon ring is 432621513.


Using the numbers 1 to 10, and depending on arrangements, it is possible to form 16- and 17-digit strings. What is the maximum 16-digit string for a “magic5-gon ring?


In attempting this problem, I choose to label the 5 inner nodes as i, j, k, l, and m.
α, β, γ, δ, and θ being the corresponding outer nodes.

Let x be the sum total of each triplet line, i.e.,

x = α + i + j = β + j + k = γ + k + l = δ + l + m = θ + m + i


First Observation:
For the string to be 16-digits, 10 has to be in the outer ring, as each number in the inner ring is included in the string twice. Next, we fill the inner ring in an iterative manner.

Second Observation:
There 9 numbers to choose from for the inner ring — 1, 2, 3, 4, 5, 6, 7, 8 and 9.
5 have to be chosen. This can be done in 9C5 = 126 ways.
According to circular permutation, if there are n distinct numbers to be arranged in a circle, this can be done in (n-1)! ways, where (n-1)! = (n-1).(n-2).(n-3)…3.2.1. So 5 distinct numbers can be arranged in 4! permutations, i.e., in 24 ways around a circle, or pentagonal ring, to be more precise.
So in all, this problem can be solved in 126×24 = 3024 iterations.

Third Observation:
For every possible permutation of an inner-ring arrangement, there can be one or more values of x (triplet line-sum) that serve as a possible contenders for a “magic” string whose triplets add up to the same number, x. To ensure this, we only need that the values of α through θ of the outer ring are distinct, different from the inner ring, with the greatest of these equal to 10.
Depending on the relative positioning of the numbers in the inner ring, one can narrow the range of x-values one might have to check for each permutation. To zero-down on such a range, let’s look at an example. Shown in the figure below is a randomly chosen permutation of number in the inner ring – 7, 2, 3, 4 and 5, in that order.


So 10, 9, 8, 6 and 1 must fill the outer circle. It’s easy to notice that the 5, 7 pair is the greatest adjacent pair. So whatever x is, it has to be at least 5 + 7 + 1 = 13 (1 being the smallest number of the outer ring). Likewise,  2, 3 is the smallest adjacent pair, so whatever x is, it can’t be any more than 2 + 3+ 10 = 15 (10 being the largest number of the outer ring). This leaves us with a narrow range of x-values to check – 13, 14 and 15.

Next, we arrange the 5 triplets in clock-wise direction starting with the triplet with the smallest number in the outer ring to form a candidate string. This exercise when done for each of the 3024 permutations will shortlist a range of candidates, of which, the maximum is chosen.

That’s all there is to the problem!

Here’s the Python Code. It executes in about a tenth of a second!

from itertools import permutations
from itertools import combinations
# array of candidate solutions empty at the beginning
record = []
# choose 5 numbers for inner cells between 1 and 9; there are 9C5 combinations
# the problem ask for a 16-digit number, so 10 is not to be included in inner cells
cells = range(1,10)
inner_cells = [map(int,comb) for comb in combinations(cells,5)]
# code to calculate min and max couple in an array
def minCouple(array):
answer = array[0]+array[-1]
for i in xrange(len(array)-1):
coupleSum = array[i] + array[i+1]
if coupleSum < answer:
answer = coupleSum
return answer
def maxCouple(array):
answer = 0
for i in xrange(len(array)-1):
if i==0:
coupleSum = array[0]+ array[-1]
if coupleSum > answer:
answer = coupleSum
coupleSum = array[i]+ array[i+1]
if coupleSum > answer:
answer = coupleSum
return answer
# Algorithm
for array in inner_cells:
pivot = array[0]
perm_array = array[1:]
perms = [map(int,perm) for perm in permutations(perm_array,4)]
for perm in perms:
checkArray = perm
outerRing = [el for el in range(1,11) if el not in checkArray]
xMax = minCouple(checkArray) + max(outerRing)
xMin = maxCouple(checkArray) + min(outerRing)
if xMax >= xMin:
for x in xrange(xMin, xMax+1):
i = checkArray[0]
j = checkArray[1]
k = checkArray[2]
l = checkArray[3]
m = checkArray[4]
alpha = x-i-j
beta = x-j-k
gamma = x-k-l
delta = x-l-m
theta = x-m-i
outerCalculated = [alpha, beta, gamma, delta, theta]
if sorted(outerCalculated) == sorted(outerRing):
a = [alpha, i, j]
b = [beta, j, k]
c = [gamma, k, l]
d = [delta, l, m]
e = [theta, m, i]
min_val = min(alpha, beta, gamma, delta, theta)
if alpha == min_val:
append = a+b+c+d+e
elif beta == min_val:
append = b+c+d+e+a
elif gamma == min_val:
append = c+d+e+a+b
elif delta == min_val:
append = d+e+a+b+c
elif theta == min_val:
append = e+a+b+c+d
l = [str(i) for i in append]
s = ''.join(l)
integer_list = int(s)
print max(record)
view raw hosted with ❤ by GitHub

Ans: 6531031914842725

Collatz Conjecture — What You Need to Know

Like many of my previous posts, this post too has something to do with a Project Euler problem. Here’s a sketch of the Colatz Conjecture.

The following iterative sequence is defined for the set of positive integers:

n → n/2 (n is even)
n → 3n + 1 (n is odd)

Using the rule above and starting with 13, we generate the following sequence:

13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1

So basically, it’s just this. Take any natural number n. If n is even, divide it by 2 to get n / 2. If n is odd, multiply it by 3 and add 1 to obtain 3n + 1. Repeat the process indefinitely. The conjecture is that no matter what number you start with, you will always eventually reach 1. The property has aptly been called oneness! But perhaps oneness has its pitfalls too…

If the conjecture is false, it can only be because there is some starting number which gives rise to a sequence that does not contain 1. Such a sequence might enter a repeating cycle that excludes 1, or increase without bound. No such sequence has been found.

It can be seen that the sequence:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1. Which starting number, under one million, produces the longest chain?

NOTE: Once the chain starts the terms are allowed to go above one million.


Histogram of stopping times for the numbers 1 to 100 million. Stopping time is on the x axis, frequency on the y axis.

Approach 1 (A naïve, but straigh forward method)

# Longest Collatz Sequence under a million
# Function listing collatz sequence for a number
def collatz(n):
"function listing collatz sequence for a positive integer"
coll = []
while n != 1:
if n % 2 == 0:
n = n/2
n = 3*n + 1
return coll
longest = 0
j = 0
for i in xrange(1, 1000000):
lencoll = len(collatz(i))
if lencoll > longest:
longest = lencoll
j = i
print j
view raw hosted with ❤ by GitHub

Approach 2 (Smart, quick method that uses dynamic programming with the help of dictionaries)

collatz = {1:1}
def Collatz(n):
global collatz
if not collatz.has_key(n):
if n%2 == 0:
collatz[n] = Collatz(n/2) + 1
collatz[n] = Collatz(3*n + 1) + 1
return collatz[n]
for j in range(1000000,0,-1):
print collatz.keys()[collatz.values().index(max(collatz.values()))]
view raw hosted with ❤ by GitHub

I couldn’t help appreciate the elegance of the second algorithm. It’ll be well worth perusing if you don’t get it at one go. [Hint: It keeps track of the number of terms of a particular sequence as values assigned to keys of a Python dictionary]

Ans: 837799

Large sum — Project Euler (Problem 13)

This isn’t much of a problem really, but since I’m posting solutions to all the Project Euler problems I solve, I’ve been OCD’d into posting this one too. Besides, it illustrates the simplifying power of Python as a language?

Anyway… here’s the problem:

Work out the first ten digits of the sum of the following one-hundred 50-digit numbers:

The solution:
I first copy the problem matrix to a .txt file, in this case, euler13.txt
The solution is cake really, and I don’t care whether this was worth posting on my blog or not coz I wasted my time solving this problem anyway, and it shouldn’t have been for nothing!
# Read the problem matrix into an array in python
filename = 'euler13.txt'
with open(filename, "r") as ins:
array = []
for line in ins:
# Convert the array into an array of integers
newArray = []
for i in array:
# Sum up the array and print the first 10 numbers of the sum as a string
arraySum = sum(newArray)
print str(arraySum)[:10]
view raw hosted with ❤ by GitHub
Ans: 5537376230