Friday, March 8, 2013

R is interesting and powerful

Recently I'm participating a online course Data Analysis by Jeff Leek , which introduces some basic ideas of data analysis. But the best part of it is through this course I finally started to use R.
Sometime at work, I need to draw some graph, usually I will use excel for simple lines or bars, or pies. and use desmos to explore functions. I also used matlab to process data.
I hear about R along ago, but there no motivation to just learn a new thing. until recently, the data analysis course lead me to explore R. It's grammar is like javascript which makes it more easier for me.
I think the reason I like to use R is totally because of RStudio which makes it a rival of matlab.

Tuesday, February 26, 2013

Practice on Matrix Factorization Using Pig

  1. Background

Recently I'm thinking about optimize my recommender system. The online version is using Item-based Collaborative Filtering, which basically recommend other items when user is viewing an item.
The Original usage of Matrix Factorization by Netflix-prize-winning algorithm is to predict rating for not viewed items and choose the highest rating items to recommend. But it take a lot of effort to predict because you have to calculate for every single item for every single user which is O(m*n). This is not practical for large system having 10s of million users and millions of items.
My variation of using Matrix Factorization is to add it on top of Item-based Collaborative Filtering. Which means:
first, using cf to generate recommend candidates, usually 10s or 100
then, for each single user, re-rank candidates according to features learned by mf.
this usage requires features learned by mf can be retrieved in real-time.  the best to do is to minimize the size of features, which means the features learned have to be as sparse as possible. I'm thinking keep about 10 significant features for each user or item.

2. Practice

according to  Matrix Factorization Techniques for Recommender System, I choose Stochastic Gradient Descend algorithm to learn the features which is:
1. easy to implement
2. relative fast running time
3. easy to parallel for large data set

the lost function is give below
lost = sum(r_ui - q_i*p_u)^2) for (i, u) + lambda*( sum(||q_i||) for i + sum(||p_u||) for u )

Problem 1, which norm should be used as regularization ?  

According to Sparse methods for machine learning Theory and algorithms
using l1-norm rather than l2-norm used in the previous paper is better as it performs regularization as well as model selection.

define: e_iu = r_ui - q_i*p_u
take partial derivative for q_i and p_u for lost function then the update for q_i and p_u are as below:

q_i = q_i + LearningRate*(sum(e_ui*p_u) for u - lambda*sign(q_i) )
p_u = p_u + LearningRate*(sum(e_ui*q_i) for i - lambda*sign(p_u) )

Problem 2, how to initiate q_i and p_u?

Choosing a right initial point is very import for stochastic gradient descent. even for a convex problem, choose a good initial point will speed up the learning process. for non-convex problem like this, a good choice also decide the quality of local optimum.
Unfortunately, there is not general rule to choose a good initial point that make a local optimum close to global optimum, but there is still way to speed up the learning process that save up a lot time.
  • Is choosing all 0 a good choice? Take a look at the update formula above, if the initial is zero, the square error part would become zero, and the regularization part will become 1. after the update, all elements of q_i and p_u will be -LearningRate*lambda. if that is what   we want, we can just initiate q_i and p_u with -LearningRate*lambda and save all calculation.
  • Is initiate every elements with a small value like -LearningRate*lambda a good choice? I don't think so, if q_i and p_u is very small, then the update of them will also be very small for each round of training, then it will take several more step to came to a point close to optimum.
My intuitive is choose the initial so that q_i*p_u is close to mean(r_ui) so that the initial is close to optimum.
Assume mean(r_ui) is 4, and q_i and p_u both are vectors with 100 elements. E() is for Expectation in statistics
E(q_i*p_u) = 100*E(q_i_k) * E(p_u_k) as q_i_k and p_u_k are independent.
E(q_i*p_u) = 100*E(q_i_k) ^2 as assuming q_i_k and p_u_k have the same distribution.
then  E(q_i_k)  = 0.2 
The question left is how to assign every q_i_k so that E(q_i_k) is 0.2.
Constant value obviously is not a good choice, because in this case, every element in a vector will update with the same difference and keep all elements the same forever.
A simple alternative is randomly pick from [0, 0.4],  or use other distribution like Gaussian Distribution with mean = 0.2. I didn't try but blindly choose random assign.

Problem 3, how to choose Learning Rate ?

Theoretically as long as the Learning Rate is small enough, stochastic gradient descent will eventually get to local optimum. so a brute force way to do it is to start from 0.1 and exponentially shrink like 0.01, 0.001... until every step lost function result decreases.
In my experiment, I have to choose 10E-5 to make lost function result decrease. and with this small learning rate, the decrease of lost function result is very very slow, after run about 100 round, it only decrease about 10%.
After think about this for a few days. I realized that the algorithm only take account of observed data, which lead to huge unbalance. As the hotness of items coincides with Zipf's Law, there is huge unbalance between items or users.
For example, item A may be viewed by 1000 users, but item B may only be viewed by 1 user.
then the square error part of update function of q_i for A is supposed to be 1000 times larger then B. In this case, a uniform learning rate for both item A and B would be inefficient. That is also why only such little learning rate can decrease lost function result. It's because after sum errors from all users view a popular item, the update difference could be very large.
To tackle the problem, I use avg instead of sum so that the difference of update would be the same size as q_i or p_u. Then the update function becomes to:

q_i = q_i + LearningRate*(avg(e_ui*p_u) for u - lambda*sign(q_i) )
p_u = p_u + LearningRate*(avg(e_ui*q_i) for i - lambda*sign(p_u) )

In this case, I only need a 0.01 learning rate and the learning speed is good, after taken 60 round, more than half error is reduced.

Problem 4, how to make feature vector sparse?


As described in background section, I hope the feature vector of user and item to be as sparse as possible so that it will require less space to store those vectors. Further, in reality, a user usually has only several interests and an item is also belong to several categories, so spareness is necessary.
The idea is simple, to make the result spare, we just have to add a regularization element that evaluate spareness to lost function.
The first thing poped into my head is Entropy.
H=sum(-p_i*log(p_i)) 
where p is the probability of each element, which is for each element i in a vector
p_i = abs(x_i)/sum(x_i)
According to information theory, Entropy is a metrics for how many bits needed to describe the vector, the more sparse the vector is the less entropy is has. then after adding this into lost function, it becomes:
lost = sum(r_ui - q_i*p_u)^2) for (i, u) + lambda_1*( sum(||q_i||) for i + sum(||p_u||) for u ) + lambda_2*( sum(H(q_i)) for i + sum(H(p_j)) for j )
let's change update function next.
to simplify, we just assume sum(x_i) is a constant then it can be ignored because magnitude is not important in gradient descend.
then
d H/d p_i = -( log(p_i) + 1 )
d H/ d x_i = -( log(x_i) + 1 ) if x_i >=0
                = log(-x_i) + 1 if x_i <0
the derivative is colored in red and orange in the picture below:
Derivative of Entropy
Let's take a look of the derivative. It becomes extreme large when the absolute value of x_i is near zero. If we apply this into gradient descend, we will subtract(gradient descend use the negative of derivative) a big number when x_i is just a little bigger then 0, which make the next x_i a negative number with a very big absolute value. This will make x_i very unstable.
Let's think about how entropy makes the vector sparse when applying gradient descend. When x_i is bigger than 0.1 which makes H(x_i)=0, it add a little bit to x_i, and the bigger x_i is, the bigger the addition is. When x_i is smaller than 0.1, it subtract a little bit to x_i, and the smaller x_i is, the bigger the subtraction is. Then we can find a function whose derivative has similar attributes and the derivative of it won't be extreme large when the absolute value is small.
The simplest form is
H(p_i) = -x^2+x 
d H/d x_i = -2x+1 if x>=0
               = -2x-1 if x<0
whose derivative is simply a linear function as showed in blue and black in the above picture.
In this case, the update functions become:

q_i = q_i + LearningRate*(avg(e_ui*p_u) for u - lambda_1*sign(q_i) - lambda_2*(-2x+sign(q_i)) )
p_u = p_u + LearningRate*(avg(e_ui*q_i) for i - lambda*sign(p_u)  - lambda_2*(-2x+sign(p_j)) )

In my experiment, the entropy really make the output vector more sparse. In the picture below is the same vector trained with and without entropy regularization. Elements are displayed from small to big. As you can see the blue one is more steeper.



Thursday, January 10, 2013

a simple Maximal Clique algorithm for map reduce

Suppose there is a 6-node-graph as the pic below,  and I will explain how my algorithm finds all Maximal Cliques from it in a map-reduce friendly way.

Section 1. find all cliques level by level

Step 0. generate all node pairs for each edge and sort node in order. for this example:
(1,2) (1,3) (1,4) (2,3) (2,4) (2,5) (2,6) (3,4) (3,5) (3,6) (4,5) (4,6) (5, 6)
actually, these are all cliques with 2 nodes.
Step 1. from here we can suppose we have all cliques with N nodes, name it CliqueNs. Then cluster each clique according to the first N-1 node. for this example:
(1: {(1,2) (1,3) (1,4)}) 
(2: {(2,3) (2,4) (2,5) (2,6) }) 
(3: (3,4) (3,5) (3,6)) 
(4: (4,5) (4,6)) (5: {(5, 6)})
Step 2. for each cluster generate all pairs of clique-n as a clique-n+1 candidate,
for this example: 
(1,2,3) (1,2,4) (1,3,4)    
(2,3,4) (2,3,5) (2,3,6) (2,4,5) (2,4,6) (2,5,6)    
(3,4,5) (3,4,6) (3,5,6)    
(4,5,6)
Step 3. map all candidate using last n-1 nodes as key, map all nodes in CliqueN using all nodes as key, and reduce candidates and CliqueN together. for this example:
((1,2): {} {(1,2)}) 
((1,3): {} {(1,3)}) 
((1,4): {} {(1,4)}) 
((2,3): {(1,2,3)} {(2,3)}) 
((2,4): {(1,2,4)} {(2,4)}) 
((2,5): {} {(2,5)}) 
((2,6): {} {(2,6)}) 
((3,4): {(1,3,4) (2,3,4)} {(3,4)}) 
((3,5): {(2,3,5)} {(3,5)}) 
((3,6): {(2,3,6)} {(3,6)}) 
((4,5): {(2,4,5) (3,4,5)} {(4,5)}) 
((4,6): {(2,4,6) (3,4,6)} {(4,6)}) 
((5,6): {(2,5,6) (3,5,6) (4,5,6)} {(5,6)}) 
Step 4. only keep the entries that both candidate set and CliqueN set is not empty, and all candidates from remained entry becomes CliqueN+1, for this example:
(1,2,3) (1,2,4) (1,3,4) (2,3,4) (2,3,5) (2,3,6) (2,4,5) (3,4,5) (2,4,6) (3,4,6) (2,5,6) (3,5,6) (4,5,6)
Step 5. repeat step 2, 3, 4 until there is no CliqueN+1 or numbers in CliqueN+1 <= n+1

For this example:
redo step 1:
((1,2): {(1,2,3) (1,2,4)})
((1,3): {(1,3,4)})
((2,3): {(2,3,4) (2,3,5) (2,3,6)})
((2,4): {((2,4,5) (2,4,6))})
((2,5): {(2,5,6)})
((3,4): {(3,4,5) (3,4,6)})
((3,5): {(3,5,6)})
((4,5): {(4,5,6)})
redo step 2:
(1,2,3,4) (2,3,4,5) (2,3,4,6) (2,3,5,6) (2,4,5,6) (3,4,5,6)
redo step 3:
((1,2,3): {} {(1,2,3)})
((1,2,4): {} {(1,2,4)})
((1,3,4): {} {(1,3,4)})
((2,3,4): {(1,2,3,4)} {(2,3,4)})
((2,3,5): {} {(2,3,5)})
((2,3,6): {} {(2,3,6)})
((2,4,5): {} {(2,4,5)})
((2,4,6): {} {(2,4,6)})
((2,5,6): {} {(2,5,6)})
((3,4,5): {(2,3,4,5)} {(3,4,5)})
((3,4,6): {(2,3,4,6)} {(3,4,6)})
((3,5,6): {(2,3,5,6)} {(3,5,6)})
((4,5,6): {(2,4,5,6) (3,4,5,6)} {(4,5,6)})
redo step 4:
(1,2,3,4) (2,3,4,5) (2,3,4,6) (2,3,5,6) (2,4,5,6) (3,4,5,6)
these are clique with 4 nodes
go on like this you will eventually find all cliques.

How to implement it in map-reduce?

As you can see, almost all operations in the steps are mapping into key-values and aggregate by keys, which is the nature of map-reduce.
Step 1 and step 2 can be done within 1 map-reduce, and step 3 and step 4 can be done within 1 map-reduce.
Then we can implement this with only 2 map-reduce for each level of cliques.

Why does this algorithm work?

It bases on the fact that if 2 n-node-clique share n-1 nodes and if the 2 not shared nodes are connected, then the union of the 2 clique is also a clique, which is a n+1 node clique.



Section 2. find all maximal cliques level by level

As you can see, above steps only finds all clique but not maximal clique. 
So what is maximal clique? 
It's a clique that is not a subset of any other clique.
Think otherwise, if a clique of N nodes is not maximal, then is must be a subset of some N+1 node clique.
Then the job is easy, enumerate all n-node-subsets of n+1-node cliques and  use them to validate n-node cliques. It's easy to implement within 1 map-reduce job as step 3 and step 4 above.


Section 3. advances and pitfalls 

Advances

1. most heavy duty is implemented by hashing like step 3, which is scalable by using map-reduce
2. the algorithm is intutive and easy to implement
3. maximal cliques are found level by level, you can skip small maximal cliques if you only cares about large ones. but it's pity that find cliques can not skip.

pitfalls

1. in step 2, which is generating candidates, will generate n^2/2 (n is the size of a set of cliques) emit for each set of cliques, which is gigantic when the graph is dense. especially when implemented using map-reduce, the reduce size is estimate by map size, but in this case, the reduce size will far larger then map size, in practice it will need manual settings of reduce size. However, real-life large scale graphs or networks are usually sparse, like social network, people only have hundreds or thousands of friends. 
2. cliques is calculated level by level, it will need m (m is the size of the largest clique) iterations to complete. Map-reduce have overhead to prepare and start, suppose 1 minutes per iteration, it will take half an hour overhead for a 30 iterations, which is not rare.

about why I'm thinking about this algorithm, see my another blog  

Tuesday, December 18, 2012

an openssl sha512 example


/*--------cpp code-------------*/

#include <openssl/sha.h>
#include <stdio.h>
int main(){
        unsigned char *md = new unsigned char[64];
        SHA512((unsigned char*)"a", 1, md);
        for(int i=0; i<64; ++i){
                printf("%02x", md[i]);
        }
        printf("\n");
}

compile using:
>g++ -I /usr/local/ssl/include -L /usr/local/ssl/lib -l crypto sha512test.c
>./a.out
1f40fc92da241694750979ee6cf582f2d5d7d28e18335de05abc54d0560e0f5302860c652bf08d560252aa5e74210546f369fbbbce8c12cfc7957b2652fe9a75

Thursday, December 6, 2012

study in community detection in network

Recently I'm thinking if I can find out any interest preference within a bunch of online novel reading logs. The idea is intuitive, people may be interested in some kind of novel like love story, magic, sci-fi or military etc. , then the the novels in the same kind may be read by a same person frequently than two novels out of the same kind. It seems easy to identify for coarse categories as I already have category tag for each of the novel, but identify interest within more subdivided categories or even within totally another kind of hierarchy of categories like interest of different ages seems to have much more fun.
The first thing comes into my head is using clustering, classical algorithms are hierarchical, k-means, canopy, Gaussian mixture model, Dirichlet process clustering, etc. The limitation is using only novel reading logs, novel can not be represented as point in space. In this case, k-means, center-based hierarchical and Gaussian mixture model is not suitable. In spite of this, normal cluster algorithm only considers links between 2 nodes, but rarely consider links among all nodes within a cluster.
Clustering Coefficient is a good way to describe the attribute of a cluster, but exhaustively enumerate all possible clusters and calculate clustering coefficient is clearly unacceptable.
the paper Finding and evaluating community structure in networks, MEJ Newman, M Girvan - Physical review E, 2004 - APS is trying to solve similar problem, and the features of its method is, first it's an dividing method not an agglomerating method, second, it uses "betweenness" which is defined on an edge as weighted sum of shorted path between any 2 nodes passes through this edge, to divide the graph.
but it's different from what I thought.
first, In real world, a node can belong to different communities at the same time
second, it's fine some nodes are alone, or the cluster is very small, but the middle size clusters which is very cohesive is what I'm looking for. so I'm still prefer aggregating method of clustering.

Recently, I'm thinking through the clique method, it is:

  • first generate all maximal cliques of the whole graph. 
  • second, use the maximal cliques as core and attach other nodes to a clique to form a cluster if the node is close enough to most of the nodes in the clique. 

like the picture below, the red nodes and deep green nodes are 2 group of maximal cliques, and the orange nodes are attached to red clique to form a cluster, and grass green nodes are attached to deep green nodes to form another cluster. there are some white nodes which are not in any cluster, and the middle node have both orange and grass green belongs to both cluster
but the problem is finding clique is NP-Complete problem, when the graph is as large as social network, sequential  method will be too slow. there are some useful thing I found:

  1. "On Computing All Maximal Cliques Distributedly",Fábio Protti, Felipe M. G. França, and Jayme Luiz Szwarcfiter , a divide and conquer method, recursively divide graphs into v1, v2, and computer cliques in v1, v2 and between v1, v2
  2. xrime, a graph algorithm package, which include a maximal clique method. It generates the neighbors of neighbors for each node, and find maximal clique in each set of neighbors of neighbors. the method is clean and tricky but to emit all adjacent list to all neighbors seems to be too heavy.
here is a map-reduce friendly algorithm I work out to find maximal clique level by level, see the blog a simple Maximal Clique algorithm for map reduce for detail



Thursday, November 22, 2012

Apache Pig, a better way to map reduce

中文版:Apache Pig, 便捷的map reduce之道
Have used Apache Pig for about a year before I start to write anything about it. 
First, you won't want to use pig or map reduce if your data is in the scale of MB or even GB, it will no more efficient than python. But for big data, it's beyond the ability of single machine, so you have to think about parallel process, which is map reduce used for, and pig is the tool that make MR more easier and efficient.  
Second, I have to say it's very handy and elegant that you can make simple statistical or data transformation with only several lines of code. eg. you want to count how many query view for each query through a search log, the code looks like below:
---pig example for qv---
log = load 'xxx' as (query);
grp = group log by query;
qv = for each grp generate group as query, COUNT(log) as qv;
------
If you write map reduce, you have to follow the routines and copy-paste tens of lines of code before you start to write your logic. eg. in java on hadoop is append at the end, because it's really too long.
Third, pig like any other script languages, is very convenient to test or run, just modify the script and run, but with java, you have to modify, compile, package src as jar and then execute.
Fourth, when using java or even streaming on hadoop, you have to run you code on real cluster to test, but with pig, you can use "-x local" to run locally with local file system, which is not only faster but also more easier to check out results than using hdfs. you can use this feature to process small data if you're willing to.
Fifth, grammar check with "-c" makes me more comfortable while coding pig, you don't have to be worried about yourself missing one letter will cause you test again or even miss to tell the bug because of lack of unit test like with python.
Sixth, pig has a very extendable API, with which you can implement your own UDF(user defined function) to deal with more complex data element process. eg. make string transform, calculate variance of a bag of data, or read/store data into/from a DIY format. further, java functions from jdk can be directly used without additional code.
---mr example for qv---
import xxxx;
...
import xxxx;
public class QueryCount extends Configured implements Tool{
public static class Map extends Mapper<LongWritable, Text, Text, IntWritable> {

@Override
public void map(LongWritable key, Text value, Context context) throws IOException, InterruptedException{

//map logic
}
}
public static class Reduce extends Reducer<Text, IntWritable, Text, IntWritable>{
@Override
public void reduce(Text key, Iterable<IntWritable> values, Context context) throws IOException, InterruptedException{
//reduce logic
}
}

@Override
public int run(String[] args) throws Exception {
Job job = new Job(getConf());
//some settings on job

boolean success = job.waitForCompletion(true);
return success?0:1;

}

public static void main(String[] args) throws Exception{
int ret = ToolRunner.run(new QueryCount(), args);
System.exit(ret);
}


}
------