& & & & & & & &本博客所有文章分类的总目录：&
开源Math.NET基础数学类库使用总目录：
　　本文开始一一介绍Math.NET的几个主要子项目的相关功能的使用。今天先要介绍的是最基本Math.NET Numerics的最基本矩阵与向量计算。
　　如果本文章资源下载不了，或者文章显示有问题，请参考&：&
1.创建Numerics矩阵与向量
　　矩阵与向量计算是数学计算的核心，因此也是Math.NET Numerics的核心和基础。
　　Math.NET包括对向量(Vector)和矩阵(Matrix)的支持，类型也很多。其主要注意点有：索引是从0开始，不支持空的向量和矩阵，也就是说维数或者长度最少为1。它也支持稀疏矩阵和非稀疏矩阵的向量类型。其矩阵有3种类型：稀疏，非稀疏，对角。这2个类在MathNet.Numerics.LinearAlgebra命名空间。由于一些技术和表示的原因，每一种数据类型都有一个实现，例如MathNet.Numerics.LinearAlgebra.Double有一个DenseMatrix类型，Matrix&T& 是抽象类型， 要通过其他方法去初始化。可以看看源码中的定义：
1 public abstract partial class Vector&T& :IFormattable, IEquatable&Vector&T&&, IList, IList&T&
where T : struct, IEquatable&T&, IFormattable
3 public abstract partial class Matrix&T& :IFormattable, IEquatable&Matrix&T&&
where T : struct, IEquatable&T&, IFormattable
&创建也很简单，可以大概看看下面这段代码，构造函数还有更多的用法，不一一演示，要自己研究下源代码，记得要引用MathNet.Numerics.LinearAlgebra命名空间：
1 //初始化一个矩阵和向量的构建对象
2 var mb = Matrix&double&.B
3 var vb = Vector&double&.B
5 //获取随机矩阵，也可以设置随机数所属的分布
6 var randomMatrix = mb.Random(2,3);
7 //向量相当于是一个一维数组,只有长度
8 var vector0 = vb.Random(3);//也可以选择分布
10 //矩阵还可以这样初始化
11 var matrix1 = mb.Dense(2,2,0.55);
12 //使用函数初始化
13 var matrix2 = mb.Dense(2,3,(i,j)=&3*i + j );
15 //对角矩阵
16 var diagMaxtrix = mb.DenseDiagonal(3,3,5);
18 Console.WriteLine("randomMatrix: "+randomMatrix.ToString());
19 Console.WriteLine("vector0: "+vector0.ToString());
20 Console.WriteLine("matrix1: "+matrix1.ToString());
21 Console.WriteLine("matrix2: "+matrix2.ToString());
22 Console.WriteLine("diagMaxtrix: "+diagMaxtrix.ToString());
24 //当然也可以直接从数组中创建
25 double[,] x = {{ 1.0, 2.0 },{ 3.0, 4.0 }};
26 var fromArray = mb.DenseOfArray(x);
28 Console.WriteLine("fromArray: "+fromArray.ToString());
结果如下，顺便说一下，Matrix和Vector对象已经对ToString进行了重载，以比较标准化的格式化字符串输出，很方便显示和观察：
1 randomMatrix: DenseMatrix 2x3-Double
2 0.785955
3 0.878987
5 vector0: DenseVector 3-Double
8 -0.182919
10 matrix1: DenseMatrix 2x2-Double
14 matrix2: DenseMatrix 2x3-Double
18 diagMaxtrix: DenseMatrix 3x3-Double
23 fromArray: DenseMatrix 2x2-Double
2.矩阵与向量的算术运算
　　Matrix和Vector都支持常见的操作运算符号：+ ，- ， * ，/ ，%等。我们可以从源码中看到部分这样的结构，限于篇幅，只简单列举几个重载操作符的方法，详细的源码在Matrix.Operators.cs文件：
1 public static Matrix&T& operator +(Matrix&T& rightSide)
return rightSide.Clone();
5 public static Matrix&T& operator -(Matrix&T& rightSide)
return rightSide.Negate();
9 public static Matrix&T& operator *(Matrix&T& leftSide, T rightSide)
return leftSide.Multiply(rightSide);
13 public static Matrix&T& operator /(T dividend, Matrix&T& divisor)
return divisor.DivideByThis(dividend);
&矩阵的相关操作是线性代数的核心和基础，而Matrix的基础功能也是非常强大的，我们看看Matrix的关于矩阵操作的相关代码，不仅包括常见矩阵分解算法，如LU，QR，Cholesky等，而且还包括一些线性方程的求解，都是可以直接通过实例方法进行的，看看抽象类的方法原型，具体的代码在Matrix.Solve.cs文件中：
1 public abstract Cholesky&T& Cholesky();
2 public abstract LU&T& LU();
3 public abstract QR&T& QR(QRMethod method = QRMethod.Thin);
4 public abstract GramSchmidt&T& GramSchmidt();
5 public abstract Svd&T& Svd(bool computeVectors = true);
6 public abstract Evd&T& Evd(Symmetricity symmetricity = Symmetricity.Unknown);
7 public void Solve(Vector&T& input, Vector&T& result)
if (ColumnCount == RowCount)
LU().Solve(input, result);
QR().Solve(input, result);
16 public void Solve(Matrix&T& input, Matrix&T& result)
if (ColumnCount == RowCount)
LU().Solve(input, result);
QR().Solve(input, result);
26 public Matrix&T& Solve(Matrix&T& input)
var x = Build.SameAs(this, ColumnCount, input.ColumnCount);
Solve(input, x);
32 public Vector&T& Solve(Vector&T& input)
var x = Vector&T&.Build.SameAs(this, ColumnCount);
Solve(input, x);
&3.矩阵计算综合例子
　　上面的一些说明可以看到一些基本的方法情况，下面有一个实际的例子，说明基本的矩阵运算情况，当然更多高级的功能不能在一篇里面一一讲到，后续还会逐步挖掘其他使用。上代码：
2 var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
3 formatProvider.TextInfo.ListSeparator = " ";
5 //创建A，B矩阵
6 var matrixA = DenseMatrix.OfArray(new[,] { { 1.0, 2.0, 3.0 }, { 4.0, 5.0, 6.0 }, { 7.0, 8.0, 9.0 } });
7 var matrixB = DenseMatrix.OfArray(new[,] { { 1.0, 3.0, 5.0 }, { 2.0, 4.0, 6.0 }, { 3.0, 5.0, 7.0 } });
9 //矩阵与标量相乘
，使用运算符
10 var resultM = 3.0 * matrixA;
11 Console.WriteLine(@"Multiply matrix by scalar using operator *. (result = 3.0 * A)");
12 Console.WriteLine(resultM.ToString("#0.00\t", formatProvider));
13 Console.WriteLine();
15 //使用Multiply相乘，结果和上面一样
16 resultM = (DenseMatrix)matrixA.Multiply(3.0);
18 //矩阵与向量相乘 右乘
19 var vector = new DenseVector(new[] { 1.0, 2.0, 3.0 });
20 var resultV = matrixA *
23 //矩阵与向量相乘 左乘 也可以使用LeftMultiply
24 resultV = vector * matrixA;
26 //2个矩阵相乘，要注意矩阵乘法的维数要求
27 resultM = matrixA * matrixB;//也可以使用Multiply方法
28 Console.WriteLine(@"Multiply matrix by matrix using operator *. (result = A * B)");
29 Console.WriteLine(resultM.ToString("#0.00\t", formatProvider));
30 Console.WriteLine();
32 //矩阵加法 使用 + ,或者Add方法
33 resultM = matrixA + matrixB;
34 resultM = (DenseMatrix)matrixA.Add(matrixB);
36 //矩阵减法 使用 - ,或者Subtract方法
37 resultM = matrixA - matrixB;
38 resultM = (DenseMatrix)matrixA.Subtract(matrixB);
40 //矩阵除法，使用 Divide
41 resultM = (DenseMatrix)matrixA.Divide(3.0);
过程比较简单，结果这里只列出部分：
1 Multiply matrix by scalar using operator *. (result = 3.0 * A)
2 DenseMatrix 3x3-Double
8 Multiply matrix by matrix using operator *. (result = A * B)
9 DenseMatrix 3x3-Double
　　资源大家可以去本系列文章的首页进行下载：
　　如果本文章资源或者显示有问题，请参考：
阅读(...) 评论()如何实现Breeze库中的DenseMatrix[BigDecimal]类型矩阵的基本运算？ - 知乎2被浏览279分享邀请回答1添加评论分享收藏感谢收起阅读(3379)
欢迎转载，转载请注明出处，徽沪一郎。
本文简要描述线性回归算法在Spark MLLib中的具体实现，涉及线性回归算法本身及线性回归并行处理的理论基础，然后对代码实现部分进行走读。
线性回归模型
机器学习算法是的主要目的是找到最能够对数据做出合理解释的模型，这个模型是假设函数，一步步的推导基本遵循这样的思路
为了找到最好的假设函数，需要找到合理的评估标准，一般来说使用损失函数来做为评估标准
根据损失函数推出目标函数
现在问题转换成为如何找到目标函数的最优解，也就是目标函数的最优化
具体到线性回归来说，上述就转换为
梯度下降法
那么如何求得损失函数的最优解，针对最小二乘法来说可以使用梯度下降法。
随机梯度下降
&如何解决这些问题呢？可以采用收缩方法(shrinkage method)，收缩方法又称为正则化(regularization)。
主要是岭回归(ridge regression)和lasso回归。通过对最小二乘估计加
入罚约束,使某些系数的估计为0。
线性回归的代码实现
上面讲述了一些数学基础，在将这些数学理论用代码来实现的时候，最主要的是把握住相应的假设函数和最优化算法是什么，有没有相应的正则化规则。
对于线性回归，这些都已经明确，分别为
Y = A*X + B 假设函数
随机梯度下降法
岭回归或Lasso法，或什么都没有
那么Spark mllib针对线性回归的代码实现也是依据该步骤来组织的代码，其类图如下所示
函数调用路径
train-&run,run函数的处理逻辑
利用最优化算法来求得最优解,optimizer.optimize
根据最优解创建相应的回归模型, createModel
runMiniBatchSGD是真正计算Gradient和Loss的地方
def runMiniBatchSGD(
data: RDD[(Double, Vector)],
gradient: Gradient,
updater: Updater,
stepSize: Double,
numIterations: Int,
regParam: Double,
miniBatchFraction: Double,
initialWeights: Vector): (Vector, Array[Double]) = {
val stochasticLossHistory = new ArrayBuffer[Double](numIterations)
val numExamples = data.count()
val miniBatchSize = numExamples * miniBatchFraction
// if no data, return initial weights to avoid NaNs
if (numExamples == 0) {
logInfo("GradientDescent.runMiniBatchSGD returning initial weights, no data found")
return (initialWeights, stochasticLossHistory.toArray)
// Initialize weights as a column vector
var weights = Vectors.dense(initialWeights.toArray)
val n = weights.size
* For the first iteration, the regVal will be initialized as sum of weight squares
* if it's L2 for L1 updater, the same logic is followed.
var regVal = pute(
weights, Vectors.dense(new Array[Double](weights.size)), 0, 1, regParam)._2
(c, v) match { case ((grad, loss), (label, features)) =&
val l = pute(features, label, bcWeights.value, Vectors.fromBreeze(grad))
(grad, loss + l)
combOp = (c1, c2) =& (c1, c2) match { case ((grad1, loss1), (grad2, loss2)) =&
(grad1 += grad2, loss1 + loss2)
* NOTE(Xinghao): lossSum is computed using the weights from the previous iteration
* and regVal is the regularization value computed in the previous iteration as well.
stochasticLossHistory.append(lossSum / miniBatchSize + regVal)
val update = pute(
weights, Vectors.fromBreeze(gradientSum / miniBatchSize), stepSize, i, regParam)
weights = update._1
regVal = update._2
logInfo("GradientDescent.runMiniBatchSGD finished. Last 10 stochastic losses %s".format(
stochasticLossHistory.takeRight(10).mkString(", ")))
(weights, stochasticLossHistory.toArray)
&上述代码中最需要引起重视的部分是aggregate函数的使用,先看下aggregate函数的定义
def aggregate[U: ClassTag](zeroValue: U)(seqOp: (U, T) =& U, combOp: (U, U) =& U): U = {
// Clone the zero value since we will also be serializing it as part of tasks
var jobResult = Utils.clone(zeroValue, sc.env.closureSerializer.newInstance())
val cleanSeqOp = sc.clean(seqOp)
val cleanCombOp = sc.clean(combOp)
val aggregatePartition = (it: Iterator[T]) =& it.aggregate(zeroValue)(cleanSeqOp, cleanCombOp)
val mergeResult = (index: Int, taskResult: U) =& jobResult = combOp(jobResult, taskResult)
sc.runJob(this, aggregatePartition, mergeResult)
aggregate函数有三个入参,一是初始值ZeroValue,二是seqOp,三为combOp.
seqOp seqOp会被并行执行,具体由各个executor上的task来完成计算
combOp combOp则是串行执行, 其中combOp操作在JobWaiter的taskSucceeded函数中被调用
为了进一步加深对aggregate函数的理解,现举一个小小例子。启动spark-shell后,运行如下代码
val z = sc. parallelize (List (1 ,2 ,3 ,4 ,5 ,6),2)
z.aggregate (0)(math.max(_, _), _ + _)
// 运 行 结 果 为 9
res0: Int = 9
仔细观察一下运行时的日志输出, aggregate提交的job由一个stage(stage0)组成,由于整个数据集被分成两个partition,所以为stage0创建了两个task并行处理。
LeastSquareGradient
讲完了aggregate函数的执行过程, 回过头来继续讲组成seqOp的pute函数。
LeastSquareGradient用来计算梯度和误差,注意cmopute中cumGraident会返回改变后的结果。这里计算公式依据的就是cost-function中的▽Q(w)
class LeastSquaresGradient extends Gradient {
override def compute(data: Vector, label: Double, weights: Vector): (Vector, Double) = {
val brzData = data.toBreeze
val brzWeights = weights.toBreeze
val diff = brzWeights.dot(brzData) - label
val loss = diff * diff
val gradient = brzData * (2.0 * diff)
(Vectors.fromBreeze(gradient), loss)
override def compute(
data: Vector,
label: Double,
weights: Vector,
cumGradient: Vector): Double = {
val brzData = data.toBreeze
val brzWeights = weights.toBreeze
//dot表示点积，是接受在实数R上的两个向量并返回一个实数标量的二元运算，它的结果是欧几里得空间的标准内积。
//两个向量的点积写作a&b。点乘的结果叫做点积，也称作数量积
val diff = brzWeights.dot(brzData) - label
//下面这句话完成y += a*x
brzAxpy(2.0 * diff, brzData, cumGradient.toBreeze)
diff * diff
在上述代码中频繁出现breeze相关的函数,你一定会很好奇,这是个什么新鲜玩艺。
说 开 了 其 实 一 点 也 不 稀 奇, 由 于 计 算 中 有 大 量 的 矩 阵(Matrix)及 向量(Vector)计算,为了更好支持和封装这些计算引入了breeze库。
Breeze, Epic及Puck是scalanlp中三大支柱性项目, 具体可参数www.scalanlp.org
正则化过程
根据本次迭代出来的梯度和误差对权重系数进行更新,这个时候就需要用上正则化规则了。也就是下述语句会触发权重系数的更新
val update = pute(
weights, Vectors.fromBreeze(gradientSum / miniBatchSize), stepSize, i, regParam)
以岭回归为例，看其更新过程的代码实现。
class SquaredL2Updater extends Updater {
override def compute(
weightsOld: Vector,
gradient: Vector,
stepSize: Double,
iter: Int,
regParam: Double): (Vector, Double) = {
// add up both updates from the gradient of the loss (= step) as well as
// the gradient of the regularizer (= regParam * weightsOld)
// w' = w - thisIterStepSize * (gradient + regParam * w)
// w' = (1 - thisIterStepSize * regParam) * w - thisIterStepSize * gradient
val thisIterStepSize = stepSize / math.sqrt(iter)
val brzWeights: BV[Double] = weightsOld.toBreeze.toDenseVector
brzWeights :*= (1.0 - thisIterStepSize * regParam)
brzAxpy(-thisIterStepSize, gradient.toBreeze, brzWeights)
val norm = brzNorm(brzWeights, 2.0)
(Vectors.fromBreeze(brzWeights), 0.5 * regParam * norm * norm)
计算出权重系数(weights)和截距intecept，就可以用来创建线性回归模型LinearRegressionModel，利用模型的predict函数来对观测值进行预测
class LinearRegressionModel (
override val weights: Vector,
override val intercept: Double)
extends GeneralizedLinearModel(weights, intercept) with RegressionModel with Serializable {
override protected def predictPoint(
dataMatrix: Vector,
weightMatrix: Vector,
intercept: Double): Double = {
weightMatrix.toBreeze.dot(dataMatrix.toBreeze) + intercept
&注意LinearRegression的构造函数需要权重(weights)和截距(intercept)作为入参，对新的变量做出预测需要调用predictPoint
一个完整的示例程序
在spark-shell中执行如下语句来亲自体验一下吧。
import org.apache.spark.mllib.regression.LinearRegressionWithSGD
import org.apache.spark.mllib.regression.LabeledPoint
import org.apache.spark.mllib.linalg.Vectors
// Load and parse the data
val data = sc.textFile("mllib/data/ridge-data/lpsa.data")
val parsedData = data.map { line =&
val parts = line.split(',')
LabeledPoint(parts(0).toDouble, Vectors.dense(parts(1).split(' ').map(_.toDouble)))
// Building the model
val numIterations = 100
val model = LinearRegressionWithSGD.train(parsedData, numIterations)
// Evaluate model on training examples and compute training error
val valuesAndPreds = parsedData.map { point =&
val prediction = model.predict(point.features)
(point.label, prediction)
val MSE = valuesAndPreds.map{case(v, p) =& math.pow((v - p), 2)}.mean()
println("training Mean Squared Error = " + MSE)
再次强调，找到对应的假设函数，用于评估的损失函数，最优化求解方法，正则化规则原文链接:
阅读排行榜3987人阅读
scala matrix breeze（2）
/scalanlp/breeze/wiki/Quickstart
/scalanlp/breeze/wiki/Linear-Algebra-Cheat-Sheet
根据(一)中讲的，我们只需要在sbt命令行下set 包以及链接，就可以引入breeze线性代数包。然后console进入scala交互行。
scala& import breeze.linalg._
scala &import breeze.numerics._
&1 创建一个全零向量
& & & scala& val x = DenseVector.zeros[Double](5)
& & & x: breeze.linalg.DenseVector[Double] = DenseVector(0.0, 0.0, 0.0, 0.0, 0.0)
&2 创建一个全一向量
scala& val x=DenseVector.ones[Double](5)
x: breeze.linalg.DenseVector[Double] = DenseVector(1.0, 1.0, 1.0, 1.0, 1.0)
3& 创建一个n个特殊值的向量
scala& val x=DenseVector.fill(6){5.0}
x: breeze.linalg.DenseVector[Double] = DenseVector(5.0, 5.0, 5.0, 5.0, 5.0, 5.0)
scala& val x=DenseVector.fill(6){5}
x: breeze.linalg.DenseVector[Int] = DenseVector(5, 5, 5, 5, 5, 5)
4 创建一个列向量
scala& val x=DenseVector(1,2,3,4)
x: breeze.linalg.DenseVector[Int] = DenseVector(1, 2, 3, 4)
5 创建一个行向量
scala& val x=DenseVector(1,2,3,4).t
x: breeze.linalg.Transpose[breeze.linalg.DenseVector[Int]] = Transpose(DenseVector(1, 2, 3, 4))
6 通过一个array创建一个向量
scala& val x=DenseVector(Array(1,2,3,4))
x: breeze.linalg.DenseVector[Int] = DenseVector(1, 2, 3, 4)
7 随机产生有n个值的向量
scala& val x=DenseVector.rand(4)
x: breeze.linalg.DenseVector[Double] = DenseVector(0.9, 0.684)
8 通过函数产生向量
scala& val x=DenseVector.tabulate(3){i=&2*i}
x: breeze.linalg.DenseVector[Int] = DenseVector(0, 2, 4)
1 产生一个全零的矩阵
scala& val m=DenseMatrix.zeros[Double](3,4)
m: breeze.linalg.DenseMatrix[Double] =
0.0& 0.0& 0.0& 0.0 &
0.0& 0.0& 0.0& 0.0 &
0.0& 0.0& 0.0& 0.0
2& 创建一个对角线全1.0的方阵
scala& val m=DenseMatrix.eye[Double](4)
m: breeze.linalg.DenseMatrix[Double] =
1.0& 0.0& 0.0& 0.0 &
0.0& 1.0& 0.0& 0.0 &
0.0& 0.0& 1.0& 0.0 &
0.0& 0.0& 0.0& 1.0 &
3 创建一个对角方阵
scala& val m=diag(DenseVector(1.0,2.0,3.0))
m: breeze.linalg.DenseMatrix[Double] =
1.0& 0.0& 0.0 &
0.0& 2.0& 0.0 &
0.0& 0.0& 3.0
4 可以这样创建矩阵，也可以这样赋值
scala& val m=DenseMatrix((1.0,2.0),(3.0,4.0))
m: breeze.linalg.DenseMatrix[Double] =
1.0& 2.0 &
5 通过函数产算矩阵
scala& val m=DenseMatrix.tabulate(3,2){case (i,j) =&i+j}
m: breeze.linalg.DenseMatrix[Int] =
6 从array创建矩阵 （有问题）
scala& val m=DenseMatrix(2,3,Array(11,12,13,21,22,23))
&console&:13: error: could not find implicit value for parameter rl: breeze.linalg.support.LiteralRow[Any,V]
&&&&&& val m=DenseMatrix(2,3,Array(11,12,13,21,22,23))
&&&&&&&&&&&&&&&&&&&&&&& ^
7 随机产生0-1直接的数 矩阵
scala& val m=DenseMatrix.rand(2,3)
m: breeze.linalg.DenseMatrix[Double] =
0.35&& 0.0448& &
0.1295&& 0.126&
&&相关文章推荐
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：19203次
排名：千里之外
(4)(1)(1)(9)Scala（1）
scala 中有一个包叫做breeze.linalg.DenseMatrix，是专门针对矩阵运算的，目前只是用到DenseMatrix这一个，暂时讲解一下他的用法
DenseMatrix即稠密矩阵与SparseMatrix（稀疏矩阵相对应）
为了写起来方便，可以这样对其重命名为：
import breeze.linalg.{DenseMatrix =& BDM}即BDM相当于DenseMatrix
一、随机生成一个矩阵
val a :BDM[Double] = BDM.rand(2,10)
val b :BDM[Double] = BDM.rand(2,10)a和b均是二维矩阵，数值是随机生成的。
二、访问矩阵中的元素
1、按行访问，这样可以获得一行的所有数据
val a1 = a(0,::)
val b1 = b(0,::)括号中第一个参数为行数，后面::代表所有的列，即这一行的所有数据
2、访问单个元素
val a5 = a.data(4)
val b7 = b.data(6)data是他的一个方法，参数为我们要访问的第几个元素
注意1、索引是从0开始的，所以第5个元素的索引值为4
注意2、矩阵的存储是按列存储的，也就是第一列之后接着是第二列的数据，并不是按照行来存储的，这一点和MATLAB一样
&&相关文章推荐
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：15577次
排名：千里之外
原创：61篇
评论：18条
(8)(8)(11)(9)(6)(5)(4)(1)(4)(2)(2)(2)(1)}