scikit-learn 回归基础 分类：机器学习Sklearn

numpy ndarray属性：shape dtype
有切片运算
当对赋值后的容器部分元素进行修改，影响原ndarray。

ndarray Array creation routines:
 可以用来构造ndarray的一些常用接口（Ones and zeros and so on）
一些ndarray的接口：
 ndarray.T(转置)
 ndarray.flat: 类似于指针，文档中将其类比于迭代器。（这当然在两种语言中可以）
  flat具有赋值功能，index与向量初始化顺序有关。
 
 ndarray.item: 行为与直接取值运算相同，结果不同在于返回python 内置类型而非array数据类型。

由于直接学习numpy比较没有针对性，现考虑利用numpy及scipy进行scikit-learn一些结果的搭建。

1.1.Generalized Linear Models
一般最小二乘：
模型的求解是求导

[python] view plain copy

print?

1. from sklearn import linear_model  
2. clf = linear_model.LinearRegression()  
3. print clf.fit([[0, 2], [1, 1], [2, 0]], [0, 1, 2])  
4. print clf.coef_  
5. print clf.intercept_  
6.   
7. print “\n” * 5  
8.   
9. import numpy as np  
10. from numpy.linalg import inv  
11. from numpy import linalg  
12.   
13. x = np.array([[0, 2], [1, 1], [2, 0]])  
14. y = np.array([0, 1, 2])  
15.   
16. def LinearRegression(X, y):  
17.  # X_bar is the matrix for [1 X]  
18.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
19.  X_bar[:,1:] = X  
20.  if linalg.matrix_rank(X_bar) == min(X_bar.shape):  
21.   X = X_bar  
22.   
23.  return np.dot(np.dot(inv(np.dot(X.T, X)), X.T), y)  
24.   
25. print LinearRegression(x, y)  

from sklearn import linear_modelclf = linear_model.LinearRegression()print clf.fit([[0, 2], [1, 1], [2, 0]], [0, 1, 2])print clf.coef_print clf.intercept_print "\n" * 5import numpy as npfrom numpy.linalg import invfrom numpy import linalgx = np.array([[0, 2], [1, 1], [2, 0]])y = np.array([0, 1, 2])def LinearRegression(X, y): # X_bar is the matrix for [1 X] X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X if linalg.matrix_rank(X_bar) == min(X_bar.shape):  X = X_bar return np.dot(np.dot(inv(np.dot(X.T, X)), X.T), y)print LinearRegression(x, y)

简单地讲，一般最小二乘要求设计阵非（渐近）线性相关。

岭回归：这部分diy代码与调用模块结果不太一样：

[python] view plain copy

print?

1. from sklearn import linear_model  
2. clf = linear_model.Ridge(alpha = 0.5)  
3. print clf.fit([[0, 0], [0, 0], [1, 1]], [0, 0.1, 1])  
4. print clf.coef_  
5. print clf.intercept_  
6.   
7. import numpy as np  
8. from numpy.linalg import inv  
9.   
10. x = np.array([[0,0], [0,0], [1,1]])  
11. y = np.array([0, 0.1, 1])  
12.   
13. def Ridge(alpha ,X, y):  
14.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
15.  X_bar[:,1:] = X  
16.  X = X_bar  
17.   
18.  return np.dot(np.dot(inv(alpha * np.eye(X.shape[0]) + np.dot(X.T, X)), X.T), y)  
19.   
20. print Ridge(0.5, x, y)  

from sklearn import linear_modelclf = linear_model.Ridge(alpha = 0.5)print clf.fit([[0, 0], [0, 0], [1, 1]], [0, 0.1, 1])print clf.coef_print clf.intercept_import numpy as npfrom numpy.linalg import invx = np.array([[0,0], [0,0], [1,1]])y = np.array([0, 0.1, 1])def Ridge(alpha ,X, y): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X X = X_bar return np.dot(np.dot(inv(alpha * np.eye(X.shape[0]) + np.dot(X.T, X)), X.T), y)print Ridge(0.5, x, y)

利用求导方法将一般最小二乘推广到岭回归，可以直接得到。

由于岭回归与约束条件最小二乘等价，模型选择是一个问题。
岭回归是一个具有超参数的方法，超参数的选择可以通过交叉验证完成，即所谓GCV（ generalized Cross-Validation）
这种交叉验证有简单的形式 https://en.wikipedia.org/wiki/Cross-validation_(statistics)
Leave-p-out cross-validation 可以进行组合数的交叉验证 利用MSE

岭回归是在矩阵“可逆”与“精确性”间选择。这里默认地选择了p=1的形式。

[python] view plain copy

print?

1. from sklearn import linear_model  
2. clf = linear_model.RidgeCV(alphas = [0.01,0.1, 1.0, 10.0])  
3. print clf.fit([[0,0], [0,0], [1,1]], [0, 0.1, 1])  
4. print clf.alpha_  
5.   
6. import numpy as np  
7. from numpy.linalg import inv, norm  
8.   
9. x = np.array([[0, 0], [0, 0], [1, 1]])  
10. y = [0, 0.1, 1]  
11. alphas = [0.01,0.1, 1.0, 1]  
12.   
13. def RidgeCV(alphas, X, y):  
14.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
15.  X_bar[:,1:] = X  
16.   
17.  def MSE(alpha):  
18.   para = np.dot(np.dot(inv(np.dot(X_bar.T, X_bar) + np.eye(X_bar.shape[1]) * alpha), X_bar.T), y)  
19.   return norm(np.dot(X_bar, para) - y)  
20.   
21.  min_alpha_MSE = []  
22.  for element in alphas:  
23.   MSE_val = MSE(element)  
24.   if not min_alpha_MSE:  
25.    min_alpha_MSE = [element, MSE_val]  
26.   elif MSE_val < min_alpha_MSE[1]:  
27.    min_alpha_MSE = [element, MSE_val]  
28.   
29.  return min_alpha_MSE  
30.   
31. print RidgeCV(alphas, x, y)  

from sklearn import linear_modelclf = linear_model.RidgeCV(alphas = [0.01,0.1, 1.0, 10.0])print clf.fit([[0,0], [0,0], [1,1]], [0, 0.1, 1])print clf.alpha_import numpy as npfrom numpy.linalg import inv, normx = np.array([[0, 0], [0, 0], [1, 1]])y = [0, 0.1, 1]alphas = [0.01,0.1, 1.0, 1]def RidgeCV(alphas, X, y): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X def MSE(alpha):  para = np.dot(np.dot(inv(np.dot(X_bar.T, X_bar) + np.eye(X_bar.shape[1]) * alpha), X_bar.T), y)  return norm(np.dot(X_bar, para) - y) min_alpha_MSE = [] for element in alphas:  MSE_val = MSE(element)  if not min_alpha_MSE:   min_alpha_MSE = [element, MSE_val]  elif MSE_val < min_alpha_MSE[1]:   min_alpha_MSE = [element, MSE_val] return min_alpha_MSEprint RidgeCV(alphas, x, y)

下面要进行关于lasso回归的讨论，其单纯的是将岭回归的二范数约束变为一范数，由于尖点不可导，
应该讨论最优化问题，下面先对scipy中的优化方式进行简单介绍。

下面是一个无约束条件的最优问题的求解代码：

[python] view plain copy

print?

1. import numpy as np  
2. from scipy.optimize import minimize  
3.   
4. def rosen(x):  
5.  return sum(100*(x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)  
6.   
7. x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])  
8. res = minimize(rosen, x0, method=’nelder-mead’, \  
9.      options = {’xtol’:1e-8, ‘disp’: True})  
10.   
11. print res.x  

import numpy as npfrom scipy.optimize import minimizedef rosen(x): return sum(100*(x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])res = minimize(rosen, x0, method='nelder-mead', \     options = {'xtol':1e-8, 'disp': True})print res.x

这里nelder-mead是指使用单纯形法进行求解。（求解线性规划的算法，线性目标函数，线性不等式约束）
虽然不知道为什么要对非线性问题使用。
options中xtol为停止迭代的自变量误差，disp是否显示拟合结果，x0初值点。

下面再见一个非线性规划的代码：

[python] view plain copy

print?

1. def func(x, sign = 1.0):  
2.  return sign * (2 * x[0] * x[1] + 2 * x[0] - x[0]**2 - 2*x[1]**2)  
3.   
4. def func_deriv(x, sign = 1.0):  
5.  dfdx0 = sign * (-2 * x[0] + 2 * x[1] + 2)  
6.  dfdx1 = sign * (2 * x[0] - 4 * x[1])  
7.   
8.  return np.array([dfdx0, dfdx1])  
9.   
10. cons = ({’type’: ‘eq’,  
11.    ’fun’: lambda x:np.array([x[0]**3 - x[1]]),  
12.    ’jac’: lambda x:np.array([3.0 * (x[0]**2), -1.0])  
13.   
14.  }, {’type’: ‘ineq’,  
15.   ’fun’: lambda x:np.array([x[1] - 1]),  
16.   ’jac’: lambda x:np.array([0.0, 1.0])  
17.  })  
18.   
19. res = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \  
20.     method = ’SLSQP’, options = {‘disp’: True, })  
21.   
22. print res.x  
23.   
24. res = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \  
25.     constraints = cons ,method = ’SLSQP’, options = {‘disp’: True, })  
26.   
27. print res.x  

def func(x, sign = 1.0): return sign * (2 * x[0] * x[1] + 2 * x[0] - x[0]**2 - 2*x[1]**2)def func_deriv(x, sign = 1.0): dfdx0 = sign * (-2 * x[0] + 2 * x[1] + 2) dfdx1 = sign * (2 * x[0] - 4 * x[1]) return np.array([dfdx0, dfdx1])cons = ({'type': 'eq',   'fun': lambda x:np.array([x[0]**3 - x[1]]),   'jac': lambda x:np.array([3.0 * (x[0]**2), -1.0]) }, {'type': 'ineq',  'fun': lambda x:np.array([x[1] - 1]),  'jac': lambda x:np.array([0.0, 1.0]) })res = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \    method = 'SLSQP', options = {'disp': True, })print res.xres = minimize(func, [-1.0, 1.0], args = (-1.0,), jac = func_deriv, \    constraints = cons ,method = 'SLSQP', options = {'disp': True, })print res.x

这里使用SLSQP（细节并不明确）jac是导函数。
下面尝试对lasso回归使用此种优化方法求解。

lasso回归代码：
这里对diy代码的约束进行了较大扩展，接近真值，但是否过拟合？
后面还对岭回归进行了非线性求解模拟，仅仅是norm不同。

[python] view plain copy

print?

1. from sklearn import linear_model  
2. clf = linear_model.Lasso(alpha = 0.1)  
3. clf.fit([[0, 0], [1, 1]], [0, 1])  
4. print clf.predict([[1, 1]])  
5. print clf.predict([[0, 0]])  
6.   
7. import numpy as np  
8. from numpy.linalg import norm  
9. from scipy.optimize import minimize  
10. import functools  
11.   
12. def func(X ,y, beta):  
13.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
14.  X_bar[:,1:] = X  
15.  return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0])  
16.   
17. #set alpha to 0.1  
18. cons = ({  
19.  ’type’: ‘ineq’,  
20.  ’fun’: lambda beta:np.array([10 - norm(beta, 1)])  
21.  },)  
22.   
23. X = np.array([[0, 0], [1, 1]])  
24. y = np.array([0, 1])  
25.   
26. func_0 = functools.partial(func, X)  
27. func_1 = functools.partial(func_0, y)  
28.   
29. #func_1 has the only para beta  
30.   
31. res = minimize(func_1, [0, 0, 0], constraints = cons, method = ‘SLSQP’,\  
32.     options = {’disp’: True})  
33.   
34. beta = res.x  
35. print beta  
36. print np.dot(np.array([1, 1, 1]), beta)  
37. print np.dot(np.array([1, 0, 0]), beta)  
38.   
39. print  
40.   
41. #use this method to ridge regression  
42. cons = ({  
43.  ’type’: ‘ineq’,  
44.  ’fun’: lambda beta:np.array([10 - norm(beta)])  
45.  },)  
46.   
47. res = minimize(func_1, [0, 0, 0], constraints = cons, method = ‘SLSQP’,\  
48.     options = {’disp’: True})  
49.   
50. beta = res.x  
51. print beta  
52. print np.dot(np.array([1, 1, 1]), beta)  
53. print np.dot(np.array([1, 0, 0]), beta)  
54.   
55. print  
56.   
57. clf = linear_model.Ridge(alpha = 0.1)  
58. clf.fit(X, y)  
59. print clf.intercept_  
60. print clf.coef_  
61.   
62. print clf.intercept_ + np.dot(clf.coef_, np.array([1, 1]))  
63. print clf.intercept_ + np.dot(clf.coef_, np.array([0, 0]))  

from sklearn import linear_modelclf = linear_model.Lasso(alpha = 0.1)clf.fit([[0, 0], [1, 1]], [0, 1])print clf.predict([[1, 1]])print clf.predict([[0, 0]])import numpy as npfrom numpy.linalg import normfrom scipy.optimize import minimizeimport functoolsdef func(X ,y, beta): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0])

#set alpha to 0.1cons = ({ 'type': 'ineq', 'fun': lambda beta:np.array([10 - norm(beta, 1)]) },)X = np.array([[0, 0], [1, 1]])y = np.array([0, 1])func_0 = functools.partial(func, X)func_1 = functools.partial(func_0, y)#func_1 has the only para betares = minimize(func_1, [0, 0, 0], constraints = cons, method = 'SLSQP',\    options = {'disp': True})beta = res.xprint betaprint np.dot(np.array([1, 1, 1]), beta)print np.dot(np.array([1, 0, 0]), beta)print#use this method to ridge regressioncons = ({ ‘type’: ‘ineq’, ‘fun’: lambda beta:np.array([10 - norm(beta)]) },)res = minimize(func_1, [0, 0, 0], constraints = cons, method = ‘SLSQP’,\    options = {‘disp’: True})beta = res.xprint betaprint np.dot(np.array([1, 1, 1]), beta)print np.dot(np.array([1, 0, 0]), beta)printclf = linear_model.Ridge(alpha = 0.1)clf.fit(X, y)print clf.intercept_print clf.coef_print clf.intercept_ + np.dot(clf.coef_, np.array([1, 1]))print clf.intercept_ + np.dot(clf.coef_, np.array([0, 0]))

对于lasso回归具有稀疏解的问题，其直观的理解为对二维等高线及约束的图像观察，大多数取约束尖点为
最优解的原因，为登高圆圆心以概率1与约束多边形中心不垂直。（http://blog.csdn.net/xidianzhimeng/article/details/20856047）
由于很可能通过lasso回归得到稀疏解，故其合适对高维数据进行特征选取。（高维因子模型降维）
具体细节以后展开。（L1-based feature selection）

lasso回归参数alpha称为稀疏度，其的选取类似于岭回归交叉验证，细节暂时不做讨论。（还有AIC BIC等选择方法）

将lasso与Ridge进行结合得到的Elastic Net方法。下面是实现：

[python] view plain copy

print?

1. from sklearn import linear_model  
2.   
3. clf = linear_model.ElasticNet(l1_ratio = 0.1)  
4. clf.fit([[108, 33], [44, 10001]], [861, 444])  
5. print clf.coef_  
6. print clf.intercept_  
7. print clf.predict([[108, 33]])  
8. print clf.predict([[44, 10001]])  
9.   
10. print  
11.   
12. import numpy as np  
13. from numpy.linalg import norm  
14. from scipy.optimize import minimize  
15. import functools  
16.   
17. def func(X, y, beta):  
18.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
19.  X_bar[:,1:] = X  
20.  return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0])  
21.   
22. alpha = 1.0  
23. l1_ratio = 0.1  
24. penalty = 10  
25.   
26. cons = ({  
27.   ’type’:‘ineq’,  
28.   ’fun’: lambda beta: penalty - alpha * l1_ratio * norm(beta, 1) - alpha * (1 - l1_ratio) * norm(beta) / 2  
29.  },)  
30.   
31. X = np.array([[108, 33], [44, 10001]])  
32. y = np.array([861, 444])  
33.   
34. func_0 = functools.partial(func, X)  
35. func_1 = functools.partial(func_0, y)  
36.   
37. res = minimize(func_1, [0, 8, 1], constraints = cons, method = ‘SLSQP’,\  
38.      options = {’disp’: True})  
39.   
40. print res.x  
41. print np.dot([1 ,108, 33], res.x)  
42. print np.dot([1 ,44, 10001], res.x)  

from sklearn import linear_modelclf = linear_model.ElasticNet(l1_ratio = 0.1)clf.fit([[108, 33], [44, 10001]], [861, 444])print clf.coef_print clf.intercept_print clf.predict([[108, 33]])print clf.predict([[44, 10001]])printimport numpy as npfrom numpy.linalg import normfrom scipy.optimize import minimizeimport functoolsdef func(X, y, beta): X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, beta) - y) / (2 * X.shape[0])alpha = 1.0l1_ratio = 0.1penalty = 10cons = ({  'type':'ineq',  'fun': lambda beta: penalty - alpha * l1_ratio * norm(beta, 1) - alpha * (1 - l1_ratio) * norm(beta) / 2 },)X = np.array([[108, 33], [44, 10001]])y = np.array([861, 444])func_0 = functools.partial(func, X)func_1 = functools.partial(func_0, y)res = minimize(func_1, [0, 8, 1], constraints = cons, method = 'SLSQP',\     options = {'disp': True})print res.xprint np.dot([1 ,108, 33], res.x)print np.dot([1 ,44, 10001], res.x)

这里predict值都很接近自变量，但是在参数方面不太统一。
但设定为不估计常数项后结果非常接近。

对于lasso回归应用于多元情况（因变量为矩阵），仅仅是将维数进行了
扩充，但在代码实现上要注意，最优化不等式约束仅仅能对向量进行，对于矩阵，要将其拉直为向量
形式（否则出错），下面的最优化diy代码对其进行了实践。
并将参数的差距与ols形式（无偏估计）进行了比较。

[python] view plain copy

print?

1. import numpy as np  
2. from numpy.random import normal, poisson, multivariate_normal  
3. from numpy.linalg import norm, inv  
4. from functools import partial  
5. from scipy.optimize import minimize  
6.   
7. X = np.zeros([5, 5])  
8. W = np.zeros([5, 5])  
9. for i in range(5):  
10.         X[i:] = normal(10, 3, 5)  
11.         W[i:] = poisson(10, 5)  
12.   
13. Y = np.dot(X, W) + multivariate_normal(np.ones(5), np.eye(5), 5)  
14.   
15.   
16. def func(X, Y, w):  
17.         W = np.zeros([5, 5])  
18.         for i in range(len(w)):  
19.                 first, second = divmod(i, 5)  
20.                 W[first][second] = w[i]  
21.   
22.         return norm(np.dot(X, W) - Y) / (2 * X.shape[0])  
23.   
24. penalty = 10  
25. alpha = 0.01  
26.   
27. def strain(penalty, alpha, w):  
28.         W = np.zeros([5, 5])  
29.         for i in range(len(w)):  
30.                 first, second = divmod(i, 5)  
31.                 W[first][second] = w[i]  
32.   
33.         return penalty - alpha * sum(np.diag(np.dot(W.T, W)) ** 0.5)  
34.   
35. strain_0 = partial(strain, penalty)  
36. strain_1 = partial(strain_0, alpha)  
37.   
38. cons = ({’type’:‘ineq’,  
39.         ’fun’:  strain_1},)  
40.   
41. func_0 = partial(func, X)  
42. func_1 = partial(func_0, Y)  
43.   
44. res = minimize(func_1, np.zeros(5 * 5), constraints = cons, method = ‘SLSQP’,\  
45.      options = {’disp’: True})  
46.   
47. W_bar = np.zeros([5, 5])  
48. for i in range(len(res.x)):  
49.         first, second = divmod(i, 5)  
50.         W_bar[first][second] = res.x[i]  
51.   
52. print norm(W)  
53. print norm(W - W_bar)  
54. print norm(W - np.dot(np.dot(inv(np.dot(X.T, X)), X.T), Y))  

import numpy as npfrom numpy.random import normal, poisson, multivariate_normalfrom numpy.linalg import norm, invfrom functools import partialfrom scipy.optimize import minimizeX = np.zeros([5, 5])W = np.zeros([5, 5])for i in range(5):        X[i:] = normal(10, 3, 5)        W[i:] = poisson(10, 5)Y = np.dot(X, W) + multivariate_normal(np.ones(5), np.eye(5), 5)def func(X, Y, w):        W = np.zeros([5, 5])        for i in range(len(w)):                first, second = divmod(i, 5)                W[first][second] = w[i]        return norm(np.dot(X, W) - Y) / (2 * X.shape[0])penalty = 10alpha = 0.01def strain(penalty, alpha, w):        W = np.zeros([5, 5])        for i in range(len(w)):                first, second = divmod(i, 5)                W[first][second] = w[i]        return penalty - alpha * sum(np.diag(np.dot(W.T, W)) ** 0.5)strain_0 = partial(strain, penalty)strain_1 = partial(strain_0, alpha)cons = ({'type':'ineq',        'fun':  strain_1},)func_0 = partial(func, X)func_1 = partial(func_0, Y)res = minimize(func_1, np.zeros(5 * 5), constraints = cons, method = 'SLSQP',\     options = {'disp': True})W_bar = np.zeros([5, 5])for i in range(len(res.x)):        first, second = divmod(i, 5)        W_bar[first][second] = res.x[i]print norm(W)print norm(W - W_bar)print norm(W - np.dot(np.dot(inv(np.dot(X.T, X)), X.T), Y))

从结果来看alpha较小时diy multi-lasso 收敛到ols,与理论相符。

下面实现了其的sklearn版本及DIY版本，注意在自定义代码中常数项的加法，对数据阵的加法与一元情况相同，
对参数阵比数据阵第二维度大（可以推出）

[python] view plain copy

print?

1. rom sklearn import linear_model  
2. clf = linear_model.MultiTaskLasso(alpha = 0.1)  
3. clf.fit([[0,0], [1,1], [2,2]], [[0,0], [1,1], [2,2]])  
4. print clf.coef_  
5. print clf.intercept_  
6.   
7. print “\n”  
8.   
9. import numpy as np  
10. from numpy.linalg import norm  
11. from scipy.optimize import minimize  
12. from functools import partial  
13. import copy  
14.   
15. def func(X, Y, w):  
16.  W = np.resize(w, [X.shape[1] + 1, Y.shape[1]])  
17.   
18.  X_bar = np.ones([X.shape[0], X.shape[1] + 1])  
19.  X_bar[:,1:] = X  
20.   
21.  return norm(np.dot(X_bar, W) - Y) / (2 * X.shape[0])  
22.   
23. penalty = 1  
24.   
25. def strain(penalty , first_dim, second_dim ,w):  
26.  W = np.resize(w ,[first_dim, second_dim])  
27.   
28.  return penalty - sum(np.diag(np.dot(W.T, W)) ** 0.5)  
29.   
30. X = np.array([[0,0], [1,1], [2,2]])  
31. Y = copy.deepcopy(X)  
32.   
33. strain_0 = partial(strain, penalty)  
34. strain_1 = partial(strain_0, X.shape[1] + 1)  
35. strain_2 = partial(strain_1, Y.shape[1])  
36.   
37. cons = ({’type’: ‘ineq’,  
38.   ’fun’: strain_2},)  
39.   
40. func_0 = partial(func, X)  
41. func_1 = partial(func_0, Y)  
42.   
43. res = minimize(func_1, np.zeros([(X.shape[1] + 1) * Y.shape[1]]),\  
44.    method = ’SLSQP’, constraints = cons, options = {‘disp’: True})  
45.   
46. print res.x.reshape(X.shape[1] + 1, Y.shape[1])  

rom sklearn import linear_modelclf = linear_model.MultiTaskLasso(alpha = 0.1)clf.fit([[0,0], [1,1], [2,2]], [[0,0], [1,1], [2,2]])print clf.coef_print clf.intercept_print "\n"import numpy as npfrom numpy.linalg import normfrom scipy.optimize import minimizefrom functools import partialimport copydef func(X, Y, w): W = np.resize(w, [X.shape[1] + 1, Y.shape[1]]) X_bar = np.ones([X.shape[0], X.shape[1] + 1]) X_bar[:,1:] = X return norm(np.dot(X_bar, W) - Y) / (2 * X.shape[0])penalty = 1def strain(penalty , first_dim, second_dim ,w): W = np.resize(w ,[first_dim, second_dim]) return penalty - sum(np.diag(np.dot(W.T, W)) ** 0.5)X = np.array([[0,0], [1,1], [2,2]])Y = copy.deepcopy(X)strain_0 = partial(strain, penalty)strain_1 = partial(strain_0, X.shape[1] + 1)strain_2 = partial(strain_1, Y.shape[1])cons = ({'type': 'ineq',  'fun': strain_2},)func_0 = partial(func, X)func_1 = partial(func_0, Y)res = minimize(func_1, np.zeros([(X.shape[1] + 1) * Y.shape[1]]),\   method = 'SLSQP', constraints = cons, options = {'disp': True})print res.x.reshape(X.shape[1] + 1, Y.shape[1])

利用np.resize()进行横向拉直拷贝。

更多了解请浏览：http://blog.csdn.net/sinat_30665603