Monday 2 April 2018

Importance of eigenvalues and eigenvectors

Hello there,
In this post, I am going to talk about the importance of eigenvalues and eigenvectors. For those who don't know what are eigenvalues and eigenvectors:
Suppose you have a square matrix  A of size n x n, λ be the eigenvalue and X be the eigenvector of size n x 1. Then the equation :
AXλX      ---(1)
holds true.
How eigenvalues and eigenvectors are calculated isn't the aim of this post. If you want to learn the way eigenvalues and eigenvectors are calculated, go to this link.

So coming back to the main topic of this post i.e the importance of eigenvalues and eigenvectors, I would like to give a definition of eigenvalues and eigenvectors in a different way.
The eigenvalue is a value and the Eigenvector is a vector hiding in a matrix. Both are useful to solve systems of linear equations and to find differentiation of vectors.  In the equation above X is a vector.
Suppose you have λ, as well as X with you for a particular matrix A. Can you find out the eigenvalue and eigenvector for the A^2 matrix? 
The answer is: 
The eigenvector of A^2 is same as the eigenvector of A. The eigenvalue for A^2 is  λ^2. The same goes for the nth power of the matrix

Ok, enough of formulae and definition. The question still remains the same. What is the use of it? How it works in the real world?
For answering these questions, let's dive deeper.

Let's understand what is a transformation.
A transformation is a change in the current shape of the object.  It can either be in single dimension or in multiple dimension.

1.

 2.
3.

Obviously, I am not talking about the transformation in picture number 1 😅😅😅. 
In the transformation in picture number 2, you can see that the direction of the blue arrow remains the same but the direction of the red arrow changes. The length of the blue arrow increases with some constant.
Eigenvectors are those vectors whose direction does not change even after the transformation. They only change their length when subjected to any transformation. The constant by which the length of these eigenvector changes is called as it's eigenvalue.
Similarly, in the picture number 3, you can see that the direction of the red vector remains same while that of the blue vector changes.












Thursday 8 March 2018

Faster Matrix Algebra for ATLAS

Hello Everyone,
In this and coming next posts, I am going to discuss the Gsoc 2018 project  I am applying for. To have a glimpse of the idea of the project go through the link.
I have already started working on this project and this is the work I have done till now on my Github repo. Go to this link and take an overview of the work and the details about the project, it's limitations or the current part being implemented.
I will be updating this blog from now on for this project and the issues related to this project. Any kind of suggestions and feedbacks are most welcome.


Update 1:
The current issue is with the speed of multiplication of two matrix of size 100 by 100, which takes around 317.748 milliseconds in SymMat class as compared to  24.72 milliseconds in case of Eigen:: Matrix class.
The Eigen Matrix class uses a method similar to BLASS GEMM. So the current plan is to implement something similar to this.


Update 2:
I have updated the symmetric.cpp file with using multithreading in case of multiplication of one SYmMat and One Eigen Matrix, both of size 100 by 100.
By doing so I was able to reduce the initial time complexity by 27% i.e now the time taken is 232.278 milliseconds.

Sunday 10 September 2017

Implementing your own power function

How many of you have used that pow(a,b) function in your codes? Well, I think most of you would have used it be it in C/C++/Java/Python or any other programing language which provide this inbuilt function in their libraries.
Ever wondered how that function works? or what is the time complexity in which it evaluates a^b ? What if you are asked to evaluate this on your own? What comes to your mind? The basic O(n) complexity solution? i.e multiplying a*a*a*a........ n times to compute? Can you think of a better approach having a time complexity of Olog(n)?
Yes, you can guess that if you are asked to do this in time complexity of O(log(n)) then we may end up applying the binary search. How? Let's see :

It is similar to a basic Divide and Conquer approach i.e dividing the problems into subparts and then solving it using the base case.
Suppose you have to find the value of 13^24.
 You can proceed in this way:
 13^24 = (13^12 )* (13^12)
 13^12 = (13^6)*(13^6)
 13^6 = (13^3)*(13^3)
 13^3 = (13^1) * (13^1) * 13 :-> as 3/2 = 1 (integer division )
 13^1 = (13^0)*(13^0) * 13
 What did you observe?
 if we have to compute a^n and if n is even we have to recursively compute a^(n/2) once and square it  and if n is odd we have to compute the square of  a^(n/2) and multiply it with a
 Therefore we are reducing the search space or say computing space by half which gives us the time complexity of O(log(n)).

int pow(int a, int n)
{
    if(n==0)
   {
      return 1;  //base condition
   }
   int temp = pow(a,n/2);
   if(n%2==0)
   {
       int ans = temp*temp;
      return ans;
   }
   else 
   {
      return (a*temp*temp);
   }
}



Sunday 23 July 2017

Heap in STL

  HEAP DATA STRUCTURE


#include <iostream>
#include <algorithm>
#include <vector>
using namespace std;
int main()
{
vector<int>v1;
  v1.push_back(10);
  v1.push_back(5);
  v1.push_back(8);
  v1.push_back(48);
  v1.push_back(20);
  v1.push_back(6);
  v1.push_back(1);


Making Heap 

make_heap(v1.begin(),v1.end());
This will convert this vector into a heap
cout<<v1.front()<<endl;
Sort Heap -> After this the operation the container is sorted and no longer a Heap
sort_heap(v1.begin(),v1.end());
vector<int>::iterator it ;
for(it=v1.begin();it!=v1.end();it++)
{
cout<<*it<<" ";
}
cout<<endl;
Insert new elements in Heap
push_heap() :- This function is used to insert elements into heap. The size of the heap is increased by 1. New element is placed appropriately in the heap.


 pop_heap() :- This function is used to delete the maximum element of the heap. The size of heap is decreased by 1. The heap elements are reorganized accordingly after this operation.
pop_heap() should be followed by a v1.pop_back()  otherwise it won't actually be popped. Coz pop_heap() just shifts the largest element to the end position . pop_back() reduces the size of the heap by one.
 v1.push_back(50);
     
    // using push_heap() to reorder elements
    push_heap(v1.begin(), v1.end());
     
    // Displaying the maximum element of heap
    // using front()
    cout << "The maximum element of heap after push is : ";
    cout << v1.front() << endl;
     
     // using pop_heap() to delete maximum element
    pop_heap(v1.begin(), v1.end());
    v1.pop_back();
     
    // Displaying the maximum element of heap
    // using front()
    cout << "The maximum element of heap after pop is : ";
    cout << v1.front() << endl;

}

output :
48
1 5 6 8 10 20 48
The maximum element of heap after push is : 50
The maximum element of heap after pop is : 6



is_heap() :- This function is used to check whether the container is heap or not. Generally, in most implementations, the reverse sorted container is considered as heap. Returns true if container is heap else returns false.

  is_heap_until() :- This function returns the iterator to the position till the container is the heap. Generally, in most implementations, the reverse sorted container is considered as heap.

Monday 17 July 2017

Learning Graphs

Representation :

Using adjacency list :


BFS and DFS 
// Program to print BFS traversal from a given source vertex. BFS(int s) 
// traverses vertices reachable from s.
#include<iostream>
#include <list>

using namespace std;

// This class represents a directed graph using adjacency list representation
class Graph
{
    int V;    // No. of vertices
    list<int> *adj;    // Pointer to an array containing adjacency lists

/*
adj = new list[V];
If you use simply list adj[V] system creates an array of list with size V, so your list constructor will get call V times.
new operator allocates memory and invokes the constructor.


So this statement list* adj= new list[V] allocates memory for V times the size of object of class list and invokes constructor of class list V times, and assigns the starting of the memory location to variable adj.


This kind of statement
x = new X[N];
assigns the result of the expression on the right hand side of the = to x.
So, you are dynamically allocating an array of V objects of type list, ans assigning a pointer to the first element to adj.

*/ 

    void DFSUtil(int v, bool visited[]);  // A function used by DFS
public:
    Graph(int V);  // Constructor
    void addEdge(int v, int w); // function to add an edge to graph
    void BFS(int s,int check);  // prints BFS traversal from a given source s
    void DFS(int v);    // DFS traversal of the vertices reachable from v 
};

Graph::Graph(int V)
{
    this->V = V;
    adj = new list<int>[V];
}

void Graph::addEdge(int v, int w)
{
    adj[v].push_back(w); // Add w to v’s list.
}

void Graph::BFS(int s,int  check)
{
    int ans =1;
    // Mark all the vertices as not visited
    bool *visited = new bool[V];
    for(int i = 0; i < V; i++)
        visited[i] = false;

    // Create a queue for BFS
    list<int> queue;

    // Mark the current node as visited and enqueue it
    visited[s] = true; if(s==check){ans =0;}
    queue.push_back(s);

    // 'i' will be used to get all adjacent vertices of a vertex
    list<int>::iterator i;

    while(!queue.empty())
    {
        // Dequeue a vertex from queue and print it
        s = queue.front();
         if(s==check){ans =0;}
        cout << s << " ";
        queue.pop_front();

        // Get all adjacent vertices of the dequeued vertex s
        // If a adjacent has not been visited, then mark it visited
        // and enqueue it
        for(i = adj[s].begin(); i != adj[s].end(); ++i)
        {
            if(!visited[*i])
            {
                visited[*i] = true;
                queue.push_back(*i);
            }
        }//cout<<endl;
    }cout<<endl<<ans<<endl;
}
void Graph::DFSUtil(int v, bool visited[])
{
    // Mark the current node as visited and print it
    visited[v] = true;
    cout << v << " ";

    // Recur for all the vertices adjacent to this vertex
    list<int>::iterator i;
    for (i = adj[v].begin(); i != adj[v].end(); ++i)
        if (!visited[*i])
            DFSUtil(*i, visited);
}

// DFS traversal of the vertices reachable from v. 
// It uses recursive DFSUtil()
void Graph::DFS(int v)
{
    // Mark all the vertices as not visited
    bool *visited = new bool[V];
    for (int i = 0; i < V; i++)
        visited[i] = false;

    // Call the recursive helper function to print DFS traversal
    DFSUtil(v, visited);
}


// Driver program to test methods of graph class
int main()
{

     int nodes, max_edges, origin, destin;
   // cout<<"Enter number of nodes: ";
    cin>>nodes;
    Graph g(nodes);
    max_edges = nodes * (nodes - 1);
    for (int i = 0; i < max_edges; i++)
    {
       
        cin>>origin>>destin;
        
        if((origin == -1) && (destin == -1))
            break;
         g.addEdge(origin, destin);
    }
   
    int verty,checky;
    cin>>verty>>checky;

    g.BFS(verty,checky);
    g.DFS(verty);
}


There is something called as topological sorting , so that you can check which job is scheduled first or next . you can check the order of dfs by  this :

// A C++ program to print topological sorting of a DAG
#include <iostream>
#include <stdio.h>
#include <list>
#include <stack>
using namespace std;
// Class to represent a graph
class Graph
{
    int V;    // No. of vertices'
    // Pointer to an array containing adjacency listsList
    list<int> *adj;
    // A function used by topologicalSort
    void topologicalSortUtil(int v, bool visited[], stack<int> &Stack);
public:
    Graph(int V);   // Constructor
     // function to add an edge to graph
    void addEdge(int v, int w);
    // prints a Topological Sort of the complete graph
    void topologicalSort();
};
Graph::Graph(int V)
{
    this->V = V;
    adj = new list<int>[V];
}
void Graph::addEdge(int v, int w)
{
    adj[v].push_back(w); // Add w to v’s list.
}
// A recursive function used by topologicalSort
void Graph::topologicalSortUtil(int v, bool visited[], 
                                stack<int> &Stack)
{
    // Mark the current node as visited.
    visited[v] = true;
    // Recur for all the vertices adjacent to this vertex
    list<int>::iterator i;
    for (i = adj[v].begin(); i != adj[v].end(); ++i)
        if (!visited[*i])
            topologicalSortUtil(*i, visited, Stack);
    // Push current vertex to stack which stores result
    Stack.push(v);
}
// The function to do Topological Sort. It uses recursive 
// topologicalSortUtil()
void Graph::topologicalSort()
{
    stack<int> Stack;
    // Mark all the vertices as not visited
    bool *visited = new bool[V];
    for (int i = 0; i < V; i++)
        visited[i] = false;
    // Call the recursive helper function to store Topological
    // Sort starting from all vertices one by one
    for (int i = 0; i < V; i++)
      if (visited[i] == false)
        topologicalSortUtil(i, visited, Stack);
    // Print contents of stack
    while (Stack.empty() == false)
    {
        cout << Stack.top() << " ";
        Stack.pop();
    }
}
// Driver program to test above functions
int main()
{
    // Create a graph given in the above diagram
    int V,i,j,k;
    scanf("%d",&V);
    Graph g(V+1);
    for(i=1;i<V;i++)
    {
        scanf("%d%d",&j,&k);
        g.addEdge(j,k);
    }
    cout << "Following is a Topological Sort of the given graph \n";
    g.topologicalSort();
    return 0;
}



There  is a way to traverse a rectangular maze using DFS :
function DFS(x, y, visited, n, m)
    if (x ≥ n OR y ≥ m)
        return
    if(x < 0 OR y < 0)
        return
    if(visisted[x][y] == True)
        return
    visited[x][y] = True
    DFS(x-1, y-1, visited, n, m)
    DFS(x-1, y, visited, n, m)
    DFS(x-1, y+1, visited, n, m)
    DFS(x, y-1, visited, n, m)
    DFS(x, y+1, visited, n, m)
    DFS(x+1, y-1, visited, n, m)
    DFS(x+1, y, visited, n, m)
    DFS(x+1, y+1, visited, n, m)






Difference between a graph and a tree :

In a graph there can be back edge , forward edge , cross edge , tree edge .
DFS of a graph gives you a tree. And if there are more components which are disconnected, then DFS give a tree forest.

1.In a tree  there is aunique  path from every single  vertex to other vertices. Also it is not necessary that degree of all the vertices of a tree is either 1 or 2.

2.There isn't any cycle in a tree.
 How to detect cycle in an undirected graph?
We can either use BFS or DFS. For every visited vertex ‘v’, if there is an adjacent ‘u’ such that u is already visited and u is not parent of v, then there is a cycle in graph. If we don’t find such an adjacent for any vertex, we say that there is no cycle.

3.If we take the sum of all the degrees, each edge will be counted twice. Hence, for a tree with n vertices and n – 1 edges, sum of all degrees should be 2 * (n – 1).

4.Tree has exactly n-1 edges while there is no such constraint for graph.










Sunday 2 July 2017

Babylonian algorithm of finding square root

Finding square root of a number whether its perfect square or not is simple with sqrt(n) function , but what if you have to implement it on your own?

 There are several ways for it. I will tell you about 2 methods:

1. Using Binary Search in O(log N)  time complexity

2. Babylonian Method of O(log(log N)).

1. 
   Using Binary Search

   float squareroot( float n)


{    
    // Base cases
    if (n == 0 || n == 1) 
       return n;

    // Do Binary Search for floor(sqrt(x))
     float start = 1, end = n, ans;   
    while (start <= end) 
    {       
    cout<<start<<"    "<<end<<endl; 
        float mid = (start + end) / 2;

        // If x is a perfect square

      /* It is very important to see that mid*mid can overflow, so rather than doing (mid*mid) do
     if(mid<=x/mid)      */

          if(mid==n/mid)
            return mid;

        // Since we need floor, we update answer when mid*mid is 
        // smaller than x, and move closer to sqrt(x)
         if(mid<=n/mid)
        {
            start = mid + 1;
            ans = mid;
        } 
        else // If mid*mid is greater than x
            end = mid - 1;        
    }
    return (ans+1);
}
   





2.  

float squareRoot(float n)
{
  /*We are using n itself as initial approximation
   This can definitely be improved */
  float x = n;
  float y = 1;
  float e = 0.000001; /* e decides the accuracy level*/
  while(x - y > e)
  {
    x = (x + y)/2;
    y = n/x;
  }
  return x;

}
   

Tuesday 20 June 2017

Using Sieve of Eratosthenes with lesser space complexity

This method is called as Segmented Sieve.

The idea of segmented sieve is to divide the range [0..n-1] in different segments and compute primes in all segments one by one. This algorithm first uses Simple Sieve to find primes smaller than or equal to √(n). 

Look I would really like to write about the concept of Segmented Sieve but I think the number of the viewers of my page is few, so what I am gonna do is post some links and code.

Go through this link

#include <iostream>
#include <math.h>
#include <vector>
using namespace std;
void simple_sieve(vector<int> &vec,int num)
{
int range = sqrt(num);
vector<bool>mark(range+1,true);
for(int i =2;i<=sqrt(range);i++) //checking prime numbers <= number^1/2
{
if(mark[i]==true)
{
for(int j=i*2;j<=range;j= j+i)
{
mark[j] = false;
}
}
}
for(int i=2;i<=range;i++)
{
if(mark[i]==true)
{
vec.push_back(i);
}
}

}

void  segmented_sieve(vector<int> &prime , int num)
{
simple_sieve(prime,num);
int lower_range = sqrt(num);
int limit  = lower_range;
int higher_range = 2*lower_range;


while(lower_range<num &&higher_range<=num)
{
vector<bool>mark(limit+1,true);
for(int i =0;i<prime.size();i++)
{
int k = (lower_range/prime[i])*prime[i];
if(k<lower_range){k = k+ prime[i];}
for(int j=k;j<=higher_range;j = j+prime[i])
{
mark[j-lower_range] = false;

}
}
for(int i=lower_range;i<=higher_range;i++)
{
if(mark[i-lower_range]==true)
{
prime.push_back(i);
}
}
lower_range = lower_range + limit;
higher_range = higher_range + limit;
if(higher_range>num){higher_range = num;}
}

}

int main()
{
int num = 1000;
std::vector<int> prime;
segmented_sieve(prime,num);
for (int i = 0; i < prime.size(); ++i)
{
cout<<prime[i]<<endl;
}
}


The usual time complexity of this algorithm is O(N log (log N)).
If you want to perform the Sieve of Eratosthenes in a complexity of O(n) , 
go through This link