[WIP] Multinomial newtoncg

This adds newton-cg solver to the linear_model.LogisticRegression class when multi_class="multinomial". Benchmarks on a few sparse datasets against L-BFGS solver are included below. Convergence times tended to be slightly lower for newton-cg and the scores were pretty much identical on these datasets.

Benchmark on 20NewsGroup dataset (11314 x 130107 sparse matrix):

![20news_bench](https://cloud.githubusercontent.com/assets/4610553/4179948/d35bf1d0-36dd-11e4-8bb5-b02299024725.jpeg)

Benchmark on Brown Corpus using TFIDF feature extraction (335 x 42453 sparse matrix):

![brown_bench](https://cloud.githubusercontent.com/assets/4610553/4179957/57b56b32-36de-11e4-93d9-70ba5f551238.jpeg)

Benchmark on Movie Review Corpus using TFIDF feature extraction (1340 x 39659 sparse matrix)

![moviereview_bench](https://cloud.githubusercontent.com/assets/4610553/4179803/13e0a806-36d7-11e4-9775-22426e48189a.jpeg)

Benchmark on Reuters Corpus using TFIDF feature extraction (7227 x 30916 sparse matrix)

![reuters_bench](https://cloud.githubusercontent.com/assets/4610553/4180034/40167f08-36e2-11e4-907f-90f7bbe29b7e.jpeg)


Equation (2) on page 3  of the cited reference is all that is required to get started. It gives a formula for the hessian-vector product.

![hv](https://cloud.githubusercontent.com/assets/4610553/4244364/234c724c-3a20-11e4-92fe-0ed3209d24f9.gif)

Where the arrow followed by r = 0 denotes evaluating the prior expression at r = 0. Here f is the multinomial loss function. If K is the number of classes and M is the number of features including intercept terms then v is a K*M length vector, w is the weight vector also of size K*M, and Hv is the hessian vector product of length K*M as well.

The R operator as defined by the reference is an element wise operation, we can apply this operator to each element in the gradient of f. The elements of grad(f) is

![dfdw](https://cloud.githubusercontent.com/assets/4610553/4244422/03a6f98e-3a21-11e4-8faa-830a457c9ca3.gif)

where yhat is a function of w and is calculated as

![yhat](https://cloud.githubusercontent.com/assets/4610553/4244510/dda6229e-3a22-11e4-91c2-b425355a3efa.gif)

![ak](https://cloud.githubusercontent.com/assets/4610553/4244511/e99cbbbc-3a22-11e4-9377-acf44a80fea7.gif)


and N is the number of observations. Although w is a vector I found it convenient to index it like a (K,M) matrix, the same applies for the vector v. The k,m variables above denotes this matrix index where k is the class index, and m is the feature index. Lower case y_k is our target for class k and yhat_k is the value of our estimator for class k. Since the targets and x values are not functions of w we can treat them as constants when applying the R operator. The R operator applied to the elements of grad(f) is then

![rdfdw](https://cloud.githubusercontent.com/assets/4610553/4244464/bb3fd778-3a21-11e4-8aa6-b1de790df207.gif)

The R operator applied to our estimator yhat is all that is left to compute. 

![ryhat](https://cloud.githubusercontent.com/assets/4610553/4244528/46694d24-3a23-11e4-9f8f-88f3f6a51c1d.gif)

Using the chain rule of calculus we get

![ryhat2](https://cloud.githubusercontent.com/assets/4610553/4244536/605f01f6-3a23-11e4-8d05-9178d15ea053.gif)
![ryhat3](https://cloud.githubusercontent.com/assets/4610553/4244540/73b2e984-3a23-11e4-9835-d2d88de3bddf.gif)

![ryhat4](https://cloud.githubusercontent.com/assets/4610553/4244548/aa6f02aa-3a23-11e4-9227-794db6a9c785.gif)


![ryhat5](https://cloud.githubusercontent.com/assets/4610553/4244546/9ee5bbc2-3a23-11e4-95af-ec16b152e673.gif)

![ryhat6](https://cloud.githubusercontent.com/assets/4610553/4244638/4c850a8e-3a25-11e4-8e43-459b78f6da1c.gif)

![ryhat7](https://cloud.githubusercontent.com/assets/4610553/4244708/4fe02e24-3a26-11e4-8f8a-bb516fe39169.gif)


The scalar r drops out of the equation so evaluating at r = 0 changes nothing. The variable yhat_i is a vector of all the class probability estimates (size K x 1) and x_i is the feature vector for observation i (size M x 1) and V is the matrix form of our vector v of size (K,M). 

In the hessp function the variable r_yhat computes the last line of the equations above but it computes it for all observations at once. Once r_yhat is computed all that needs to be done is a dot product with X to arrive at Hv.