最小二乘法拟合线条的C++实现

最新推荐文章于 2023-07-31 22:39:24 发布

GoodluckTian

最新推荐文章于 2023-07-31 22:39:24 发布

阅读量4.5k

点赞数 12

分类专栏：不知道放哪儿文章标签： c++ 算法机器学习

本文为博主原创文章，未经博主允许禁止转载。如有问题，欢迎指正！

本文链接：https://blog.csdn.net/weixin_43152152/article/details/119209099

版权

不知道放哪儿专栏收录该内容

5 篇文章

订阅专栏

该博客详细介绍了如何使用C++实现最小二乘法来拟合直线和二次多项式曲线，包括直接公式计算和矩阵求解方法。文章对比了不同求解方式的效率，如矩阵的SVD、QR和Cholesky分解，并提供了实际测试案例及耗时。此外，还探讨了一次和二次多项式的公式推导过程。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

最小二乘法拟合线条的C++实现

包括直线和曲线

0 参考文章

最小二乘法多项式曲线拟合数学原理及其C++实现
同事的博客，强推！原理写的很清楚！！！
Eigen官方线性算法
线性代数求解，以及矩阵求解时候的精度速度对比，可供参考
Eigen官方最小二乘对比法总结
Eigen给出的三种求解最小二乘问题的例子
最小二乘法及其c++的实现
拟合直线的公式推导及实现。

1 实现代码

总共实现了五种方法

LeastSquareFitBeeline使用推导方程式的结果公式直接计算，速度最快，精度可以调整。只能计算直线！
LeastSquareFitCurve使用推到公式直接计算二次多项式曲线。
LeastSquareNormalMethod为博客https://blog.shipengx.com/archives/e0ebe48c.html的实现过程，利用矩阵求逆计算
LeastSquareSVDMethod使用Eigen求解矩阵SVD实现
LeastSquareQRMethod使用Eigen求解矩阵的基于Householder变换的QR分解实现
LeastSquareLDLTMethod使用Eigen求解矩阵的改进型Cholesky分解实现

以上五种，第一种仅能计算直线方程，第二种计算二次多项式曲线，其余四种可以计算直线和曲线，输入相应参数即可。
正规方程LDLT、householderQr分解和bdcSvd分解，有Eigen官方给出的速度精度对比，大小数据下的速度都是不一样的
矩阵常规法是直接求逆矩阵，速度视具体数据而定
LeastSquareFitBeeline()在拟合直线时最快，因为使用的是推导的公式
选取自己最合适的方法就行了～

linefit.h

/*
 * @Author: leox_tian
 * @Date: 2021-07-30 13:06:02
 * @Description: 最小二乘法拟合线条
 */

#pragma once

#include <Eigen/Core>
#include <Eigen/Dense>
#include <chrono>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iostream>

/**参考链接：
 * http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
 * http://eigen.tuxfamily.org/dox/group__LeastSquares.html
 * https://blog.shipengx.com/archives/e0ebe48c.html
 *
 * 求解速度：正规方程LDLT>householderQr分解>矩阵常规法>bdcSvd分解
 * 求解精度：bdcSvd分解>householderQr分解>正规方程LDLT>矩阵常规法
 * LeastSquareFitBeeline()在拟合直线时最快，因为使用的是推导的公式
 */

class LineFit {
 private:
  Eigen::MatrixXf A_;
  Eigen::MatrixXf b_;

 public:
  LineFit();
  ~LineFit();
  Eigen::Vector2f LeastSquareFitBeeline(const std::vector<float>& X,
                             const std::vector<float>& Y);
  Eigen::Vector3f LeastSquareFitCurve(const std::vector<float>& X,
                                      const std::vector<float>& Y);
  void ConstructAb(const std::vector<float>& X, const std::vector<float>& Y,
                   uint8_t orders, Eigen::MatrixXf& A, Eigen::MatrixXf& b);

  Eigen::VectorXf LeastSquareNormalMethod(const std::vector<float>& X,
                                          const std::vector<float>& Y,
                                          uint8_t orders);
  Eigen::VectorXf LeastSquareSVDMethod(const std::vector<float>& X,
                                       const std::vector<float>& Y,
                                       uint8_t orders);
  Eigen::VectorXf LeastSquareQRMethod(const std::vector<float>& X,
                                      const std::vector<float>& Y,
                                      uint8_t orders);
  Eigen::VectorXf LeastSquareLDLTMethod(const std::vector<float>& X,
                                        const std::vector<float>& Y,
                                        uint8_t orders);
  Eigen::VectorXf SolveBySVD(const Eigen::MatrixXf& A,
                             const Eigen::MatrixXf& b);
  Eigen::VectorXf SolveByQR(const Eigen::MatrixXf& A, const Eigen::MatrixXf& b);
  Eigen::VectorXf SolveByLDLT(const Eigen::MatrixXf& A,
                              const Eigen::MatrixXf& b);
};

linefit.cc

#include "linefit.h"

LineFit::LineFit() {}

LineFit::~LineFit() {}

Eigen::Vector2f LineFit::LeastSquareFitBeeline(const std::vector<float>& X,
                                               const std::vector<float>& Y) {
  if (X.size() != Y.size() || X.size() < 2 || Y.size() < 2) {
    exit(EXIT_FAILURE);
  }
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

  int n = X.size();
  float sum_xy = 0.0f;
  float sum_x = 0.0f;
  float sum_y = 0.0f;
  float sum_xx = 0.0f;
  for (int i = 0; i < n; ++i) {
    sum_x += X[i];
    sum_y += Y[i];
    sum_xx += X[i] * X[i];
    sum_xy += X[i] * Y[i];
  }
  float denominator = n * sum_xx - sum_x * sum_x;
  if (fabs(denominator) <= 1e-3) denominator = 1e-3;
  // y = kx+b
  float K = (n * sum_xy - sum_x * sum_y) / denominator;
  float B = (sum_y * sum_xx - sum_x * sum_xy) / denominator;

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());

  return Eigen::Vector2f(B, K);
}
Eigen::Vector3f LineFit::LeastSquareFitCurve(const std::vector<float>& X,
                                             const std::vector<float>& Y) {
  if (X.size() != Y.size() || X.size() < 2 || Y.size() < 2) {
    exit(EXIT_FAILURE);
  }
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();
  int n = X.size();
  float sum_xy = 0.0f;
  float sum_x = 0.0f;
  float sum_y = 0.0f;
  float sum_xx = 0.0f;

  float sum_x4 = 0.0f;
  float sum_x3 = 0.0f;
  float sum_x2y = 0.0f;
  for (int i = 0; i < n; ++i) {
    sum_x += X[i];
    sum_y += Y[i];
    sum_xx += X[i] * X[i];
    sum_xy += X[i] * Y[i];
    sum_x4 += std::pow(X[i], 4);
    sum_x3 += pow(X[i], 3);
    sum_x2y += X[i] * X[i] * Y[i];
  }
  float denominator = pow(sum_xx, 3) - 2 * sum_xx * sum_x3 * sum_x -
                      sum_xx * sum_x4 * n + sum_x3 * sum_x3 * n +
                      sum_x4 * sum_x * sum_x;
  if (fabs(denominator) <= 1e-3) denominator = 1e-3;
  float a = -sum_x2y * (n * sum_xx - sum_x * sum_x) -
            sum_xy * (sum_xx * sum_x - sum_x3 * n) +
            sum_y * (sum_xx * sum_xx - sum_x3 * sum_x);
  float b = -sum_x2y * (sum_x * sum_xx - sum_x3 * n) +
            sum_xy * (sum_xx * sum_xx - sum_x4 * n) -
            sum_y * (sum_xx * sum_x3 - sum_x4 * sum_x);
  float c = sum_x2y * (sum_xx * sum_xx - sum_x3 * sum_x) -
            sum_xy * (sum_xx * sum_x3 - sum_x4 * sum_x) -
            sum_y * (sum_xx * sum_x4 - sum_x3 * sum_x3);

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());
  // y = ax²+bx+c
  c = c / denominator;
  b = b / denominator;
  a = a / denominator;
  return Eigen::Vector3f(c, b, a);
}
Eigen::VectorXf LineFit::SolveBySVD(const Eigen::MatrixXf& A,
                                    const Eigen::MatrixXf& b) {
  using namespace Eigen;
  return A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);
}
Eigen::VectorXf LineFit::SolveByQR(const Eigen::MatrixXf& A,
                                   const Eigen::MatrixXf& b) {
  return A.colPivHouseholderQr().solve(b);
}
Eigen::VectorXf LineFit::SolveByLDLT(const Eigen::MatrixXf& A,
                                     const Eigen::MatrixXf& b) {
// If the matrix A is ill-conditioned, then this is not a good method,
// because the condition number of ATA is the square of the condition number of
// A. This means that you lose twice as many digits using normal equation than
// if you use the other methods.
  return (A.transpose() * A).ldlt().solve(A.transpose() * b);
}
/**
 * @description: Ax=b，构建矩阵A和b
 * @param {XY 为对应坐标系的点坐标}
 * @param {orders 代表几次方程，为2时求解二项式曲线y=t0+t1*x+t2*x^2}
 * @return {*}
 */
void LineFit::ConstructAb(const std::vector<float>& X,
                          const std::vector<float>& Y, uint8_t orders,
                          Eigen::MatrixXf& A, Eigen::MatrixXf& b) {
  // abnormal input verification
  if (X.size() < 2 || Y.size() < 2 || X.size() != Y.size() || orders < 1)
    exit(EXIT_FAILURE);

  // map sample data from STL vector to eigen vector
  std::vector<float> x_copy(X), y_copy(Y);
  Eigen::Map<Eigen::VectorXf> sampleX(x_copy.data(), X.size());
  Eigen::Map<Eigen::VectorXf> sampleY(y_copy.data(), Y.size());

  // Vandermonde matrix of X-axis coordinate vector of sample data
  Eigen::MatrixXf mtxVandermonde(X.size(), orders + 1);

  // Vandermonde column
  Eigen::VectorXf colVandermonde = sampleX;

  // construct Vandermonde matrix column by column
  for (size_t i = 0; i < orders + 1; ++i) {
    if (0 == i) {
      mtxVandermonde.col(0) = Eigen::VectorXf::Constant(X.size(), 1, 1);
      continue;
    }
    if (1 == i) {
      mtxVandermonde.col(1) = colVandermonde;
      continue;
    }
    colVandermonde = colVandermonde.array() * sampleX.array();
    mtxVandermonde.col(i) = colVandermonde;
  }

  // calculate coefficients vector of fitted polynomial
  A = mtxVandermonde.transpose() * mtxVandermonde;
  b = mtxVandermonde.transpose() * sampleY;
}
Eigen::VectorXf LineFit::LeastSquareNormalMethod(const std::vector<float>& X,
                                                 const std::vector<float>& Y,
                                                 uint8_t orders) {
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

  ConstructAb(X, Y, orders, A_, b_);
  Eigen::VectorXf result = A_.inverse() * b_;

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());

  return result;
}

Eigen::VectorXf LineFit::LeastSquareSVDMethod(const std::vector<float>& X,
                                              const std::vector<float>& Y,
                                              uint8_t orders) {
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

  ConstructAb(X, Y, orders, A_, b_);

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());

  return SolveBySVD(A_, b_);
}
Eigen::VectorXf LineFit::LeastSquareQRMethod(const std::vector<float>& X,
                                             const std::vector<float>& Y,
                                             uint8_t orders) {
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

  ConstructAb(X, Y, orders, A_, b_);

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());

  return SolveByQR(A_, b_);
}

Eigen::VectorXf LineFit::LeastSquareLDLTMethod(const std::vector<float>& X,
                                               const std::vector<float>& Y,
                                               uint8_t orders) {
  std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

  ConstructAb(X, Y, orders, A_, b_);

  std::chrono::steady_clock::time_point t2 = std::chrono::steady_clock::now();
  std::chrono::microseconds time_used =
      std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1);
  printf("time used [%s] : %ld ms\n", __FUNCTION__, time_used.count());

  return SolveByLDLT(A_, b_);
}

2 测试demo及耗时

测试拟合直线的案例

main.cc

/*
 * @Author: leox_tian
 * @Date: 2021-07-30 13:06:02
 * @Description: mian.cc
 */
#include <algorithm>
#include <cstring>
#include <iostream>

#include "linefit.h"
using namespace std;

int main(void) {
  LineFit linefit;
  vector<float> X;
  vector<float> Y;

  cout << "一次多项式:y=5+7x" << endl;
  for (size_t i = 1; i <= 10; i++) {
    X.emplace_back(i);
    Y.emplace_back(5 + 7 * i + rand()%1) ;
  }
  int orders = 1;
  cout << linefit.LeastSquareFitBeeline(X, Y).transpose() << endl;
  cout << linefit.LeastSquareSVDMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareQRMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareLDLTMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareNormalMethod(X, Y, orders).transpose() << endl;

  cout << "\n二次多项式:y=3+x+2x^2" << endl;
  X.clear();
  Y.clear();
  for (size_t i = 1; i <= 5; i++) {
    X.emplace_back(i);
    Y.emplace_back(3 + i + 2 * i * i);
  }
  orders = 2;
  cout << linefit.LeastSquareFitCurve(X, Y).transpose() << endl;
  cout << linefit.LeastSquareSVDMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareQRMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareLDLTMethod(X, Y, orders).transpose() << endl;
  cout << linefit.LeastSquareNormalMethod(X, Y, orders).transpose() << endl;

  return 0;
}

耗时及结果
测试用例仅小数据，笔记本简单测试的，不代表实际的耗时和精度。
求解时其实很快，在构建Ab矩阵比较耗时。

一次多项式:y=5+7x
time used [LeastSquareFitBeeline] : 1 ms
5 7
time used [LeastSquareSVDMethod] : 55 ms
4.99998       7
time used [LeastSquareQRMethod] : 8 ms
5.00004 6.99999
time used [LeastSquareLDLTMethod] : 5 ms
5 7
time used [LeastSquareNormalMethod] : 19 ms
4.99998       7

二次多项式:y=3+x+2x^2
time used [LeastSquareFitCurve] : 52 ms
3 1 2
time used [LeastSquareSVDMethod] : 10 ms
2.99942  1.0004 1.99994
time used [LeastSquareQRMethod] : 7 ms
 2.9998 1.00015 1.99998
time used [LeastSquareLDLTMethod] : 4 ms
2.76033 1.18035 1.97201
time used [LeastSquareNormalMethod] : 15 ms
2.99988 1.00012 1.99998

3 关于简单公式推导

关于上面的LeastSquareFitBeeline和LeastSquareFitCurve的公式推导。
一次多项式 $y = a x + b$ 和二次多项式 $y=ax^2+bx+c$ 相对比较简单，误差平方和求偏导数后，可以直接对方程式组求解，将解的通用表达式写到代码里，也能求得最小二乘法的结果，比较笨的一种方法，不能通用求解多次多项式，但是比较简单直接。
一次多项式是可以直接手算得出通解的，二次多项式比较麻烦，所以统一用Python进行方程式组的推导计算
关于方程式组的推导计算可参考博客：Python-sympy科学计算与数据处理(方程，微分，微分方程，积分),用Python对数学函数进行求值、求偏导

以下使用juyter notebook直接导出markdown，原文件放在项目代码里，项目代码链接位于本文最后。

from sympy import *
init_printing()

a, b, c, n = symbols('a b c n')
sum_xy = symbols('\sum_{i=0}^{n}{x_iy_i}')
sum_xx = symbols('\sum_{i=0}^{n}{x_i^2}')
sum_x = symbols('\sum_{i=0}^{n}{x_i}')
sum_y = symbols('\sum_{i=0}^{n}{y_i}')
sum_x4 = symbols('\sum_{i=0}^{n}{x_i^4}')
sum_x3 = symbols('\sum_{i=0}^{n}{x_i^3}')
sum_x2y = symbols('\sum_{i=0}^{n}{x_i^2y_i}')

a, b, c, sum_xy, sum_xx, sum_x, sum_y, sum_x4,sum_x3, sum_x2y

$\left ( a, \quad b, \quad c, \quad \sum_{i=0}^{n}{x_iy_i}, \quad \sum_{i=0}^{n}{x_i^2}, \quad \sum_{i=0}^{n}{x_i}, \quad \sum_{i=0}^{n}{y_i}, \quad \sum_{i=0}^{n}{x_i^4}, \quad \sum_{i=0}^{n}{x_i^3}, \quad \sum_{i=0}^{n}{x_i^2y_i}\right )$

3.1 最小二乘法拟合直线

一次多项式 $f (x) = a x + b$
误差平方和 $s=\sum_{i=0}^{n}{(f(x_i)-y_i)^2}$
对误差平方和s求偏导,得到如下方程式组，求解ab
推导的公式简单，可以放入程序里直接使用，速度很快。

eq1 = Eq( sum_xx * a + sum_x * b - sum_xy,0)
eq2 = Eq( sum_x * a + n * b - sum_y, 0)
eq1,eq2

$\left ( \sum_{i=0}^{n}{x_i^2} a - \sum_{i=0}^{n}{x_iy_i} + \sum_{i=0}^{n}{x_i} b = 0, \quad \sum_{i=0}^{n}{x_i} a - \sum_{i=0}^{n}{y_i} + b n = 0\right )$

# aa = solve([sum_xx * a + sum_x * b - sum_xy, sum_x * a + n * b - sum_y], a,b)
aa = solve((eq1, eq2), a, b)
aa

$\left \{ a : \frac{\sum_{i=0}^{n}{x_iy_i} n - \sum_{i=0}^{n}{x_i} \sum_{i=0}^{n}{y_i}}{\sum_{i=0}^{n}{x_i^2} n - \left(\sum_{i=0}^{n}{x_i}\right)^{2}}, \quad b : \frac{\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{y_i} - \sum_{i=0}^{n}{x_iy_i} \sum_{i=0}^{n}{x_i}}{\sum_{i=0}^{n}{x_i^2} n - \left(\sum_{i=0}^{n}{x_i}\right)^{2}}\right \}$

3.2 最小二乘法拟合曲线

二次多项式 $f (x) = a x ² + b x + c$
误差平方和 $s=\sum_{i=0}^{n}{(f(x_i)-y_i)^2}$
对误差平方和s求偏导,得到如下方程式组，求解abc
到二次多项式时，abc的解已经很长了，不是很适合直接将公式放到程序里运算(当然也可以，会比较臃肿)，此时使用矩阵计算求Ax=b的方法会比较合适。

eq3 = Eq(sum_x4 * a + sum_x3 * b + sum_xx * c - sum_x2y, 0)
eq4 = Eq(sum_x3 * a + sum_xx * b + sum_x * c - sum_xy, 0)
eq5 = Eq(sum_xx * a + sum_x * b + n * c - sum_y, 0)
eq3, eq4, eq5

$\left ( - \sum_{i=0}^{n}{x_i^2y_i} + \sum_{i=0}^{n}{x_i^2} c + \sum_{i=0}^{n}{x_i^3} b + \sum_{i=0}^{n}{x_i^4} a = 0, \quad \sum_{i=0}^{n}{x_i^2} b + \sum_{i=0}^{n}{x_i^3} a - \sum_{i=0}^{n}{x_iy_i} + \sum_{i=0}^{n}{x_i} c = 0, \quad \sum_{i=0}^{n}{x_i^2} a + \sum_{i=0}^{n}{x_i} b - \sum_{i=0}^{n}{y_i} + c n = 0\right )$

bb = solve((eq3, eq4, eq5), a, b, c)
bb

$\left \{ a : \frac{- \sum_{i=0}^{n}{x_i^2y_i} \left(\sum_{i=0}^{n}{x_i^2} n - \left(\sum_{i=0}^{n}{x_i}\right)^{2}\right) - \sum_{i=0}^{n}{x_iy_i} \left(\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i} - \sum_{i=0}^{n}{x_i^3} n\right) + \sum_{i=0}^{n}{y_i} \left(\left(\sum_{i=0}^{n}{x_i^2}\right)^{2} - \sum_{i=0}^{n}{x_i^3} \sum_{i=0}^{n}{x_i}\right)}{\left(\sum_{i=0}^{n}{x_i^2}\right)^{3} - 2 \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^3} \sum_{i=0}^{n}{x_i} - \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^4} n + \left(\sum_{i=0}^{n}{x_i^3}\right)^{2} n + \sum_{i=0}^{n}{x_i^4} \left(\sum_{i=0}^{n}{x_i}\right)^{2}}, \quad b : \frac{- \sum_{i=0}^{n}{x_i^2y_i} \left(\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i} - \sum_{i=0}^{n}{x_i^3} n\right) + \sum_{i=0}^{n}{x_iy_i} \left(\left(\sum_{i=0}^{n}{x_i^2}\right)^{2} - \sum_{i=0}^{n}{x_i^4} n\right) - \sum_{i=0}^{n}{y_i} \left(\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^3} - \sum_{i=0}^{n}{x_i^4} \sum_{i=0}^{n}{x_i}\right)}{\left(\sum_{i=0}^{n}{x_i^2}\right)^{3} - 2 \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^3} \sum_{i=0}^{n}{x_i} - \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^4} n + \left(\sum_{i=0}^{n}{x_i^3}\right)^{2} n + \sum_{i=0}^{n}{x_i^4} \left(\sum_{i=0}^{n}{x_i}\right)^{2}}, \quad c : \frac{\sum_{i=0}^{n}{x_i^2y_i} \left(\left(\sum_{i=0}^{n}{x_i^2}\right)^{2} - \sum_{i=0}^{n}{x_i^3} \sum_{i=0}^{n}{x_i}\right) - \sum_{i=0}^{n}{x_iy_i} \left(\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^3} - \sum_{i=0}^{n}{x_i^4} \sum_{i=0}^{n}{x_i}\right) - \sum_{i=0}^{n}{y_i} \left(\sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^4} - \left(\sum_{i=0}^{n}{x_i^3}\right)^{2}\right)}{\left(\sum_{i=0}^{n}{x_i^2}\right)^{3} - 2 \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^3} \sum_{i=0}^{n}{x_i} - \sum_{i=0}^{n}{x_i^2} \sum_{i=0}^{n}{x_i^4} n + \left(\sum_{i=0}^{n}{x_i^3}\right)^{2} n + \sum_{i=0}^{n}{x_i^4} \left(\sum_{i=0}^{n}{x_i}\right)^{2}}\right \}$