# C++ performance: Linear regression in other way

Here is the code that can be used for calculation of mathematical function, like ax^2 + bx + c.

It is fast enough if you choose small length, otherwise if programmer don't know the small range, that code can be really slow. I've made it specially on C++ to be more fast.

#include <iostream>
#include <vector>

using namespace std;

template<class var>
var Module(var x){
if (x >= 0)
return x;
else
return x*-1;
}

class Linear {
public:
float resA, resB, resC;
float err;

float Predict(float a, float b, float c, float x) {
return ((a * (x*x)) + (b*x) + c);
}

float Predict(float x) {
return ((resA * (x * x)) + (resB * x) + resC);
}

float ErrorAv(float a, float b, float c, vector<float> input, vector<float> output) {
float error = Module(Predict(a, b, c, input[0]) - output[0]);
for (int i = 1; i < input.size(); i++)
error = (Module(Predict(a, b, c, input[i]) - output[i]) + error)/2;
return error;
}

void LinearRegr(vector<float> input, vector<float> output, float maximum, float minimum = 0, float step = 1) {
if (step == 0)
step++;
float a, b, c;

float lastError = INFINITY;
for (a = minimum; a <= maximum; a += step)
for (b = minimum; b <= maximum; b += step)
for (c = minimum; c <= maximum; c += step) {
float error = ErrorAv(a, b, c, input, output);
if (error < lastError) {
lastError = error;
resA = a;
resB = b;
resC = c;
err = lastError;

if (!lastError)
return;
}
}

}
};

#include <ctime>
int main(){
vector<float> input; //Input example.
vector<float> output; //Output example.

float a = 10.5, b = -7, c = 5.5; //Variables as search values.

//Fill dataset:
for (int i = 0; i < 100; i++) {
input.push_back(i);
output.push_back((a * (i * i)) + (b * i) + c);
}

clock_t begin = clock(); //Start clock while searching a, b, c values
Linear linear;
linear.LinearRegr(input, output, 15, -10, 0.5); //Start searching.
cout << "Time: ~" << double(clock() - begin) / CLOCKS_PER_SEC << " seconds." << endl;

cout << linear.resA << "*x^2 + " << linear.resB << "x + " << linear.resC << endl; //Enter results.
}


As can see, Linear.LinearRegr function can get from 3 to 5 parameters. I dont need make the code super prettier(for me its already pretty), just want to work it faster.

How to optimize and make it faster?

Couple suggestions

### Just use std::abs

It's in the cmath header. It's less about performance, and more about readability. Though it will likely have the same or better performance.

### Be careful with float, double, and integer constants

When working with float, make sure to append f to the end. Like 1.0f. If you use double, then you don't need to worry about this. Also, more than likely your target computer has little to no discernible performance difference between float and double

Also, be careful with = comparisons. When using float/double and doing equality operations, you run the issue of having an error and never actually getting perfect equality. Using ! on the float is fine, because you should eventually hit 0 and that will trigger the break.

// this
for (c = minimum; c <= maximum; c += step)
// should be
for (c = minimum; c < maximum; c += step)


### Use a better benchmark

I recommend Google Benchmark or quick-bench.com (it's GBench on the backend). You'll be able to see what differences your code is actually making, since it reruns the test to get a more consistent number. Google Benchmark has a lot of fun features so you can see what parameters affect your output

### Use const & for your input vectors

If you're not changing the vectors, make them const refs so the compiler knows it can take shortcuts

By rewriting ErrorAv to use std::transform_reduce when I load everything in Godbolt.org, the function appears to be inline

    float ErrorAv(float a, float b, float c, vector<float> const & input, vector<float> const & output) {
float error = abs(Predict(a, b, c, input[0]) - output[0]);
return std::transform_reduce(std::next(cbegin(input)),cend(input),std::next(cbegin(output)),error,
[](auto const & error, auto const & val){
return (val + error)/2.0;
}, // "sum function"
[this,a,b,c](auto const & i, auto const & o){
return abs(Predict(a, b, c, i) - o);
}); // "product function"
}


Since std::transform_reduce is a little strange, I'll try to explain it. Normally it takes two containers, multiplies corresponding elements together, then sums all the products together. You can overload this with lambdas, which is what I'm presenting. std::transform_reduce can be massively parallelized and take advantage of SIMD instructions, so it should provide noticeable speed up. std::transform_reduce lives in the numeric header (for good reason ;)). You'll need to compile for C++17, and if you can't, you can use std::inner_product

### Misc

• Remove the last err = lastError; in your if-statement, in the hot-path of LinearRegr. It's not doing anything.

• Put the uninitialized floats next to their loop initialization.

// this
float a,b,c;
for (a = minimum; a <= maximum; a += step)
for (b = minimum; b <= maximum; b += step)
for (c = minimum; c <= maximum; c += step)
// should be
for (float a = minimum; a <= maximum; a += step)
for (float b = minimum; b <= maximum; b += step)
for (float c = minimum; c <= maximum; c += step)


### Final note

While Rosetta code doesn't always have the best code, check out their polynomial regression algorithm which runs in linear time.

# Optimizing the algorithm

If $$\M\$$ is the number of steps between minimum and maximum, and $$\N\$$ is the size of the input and output vectors, then your algorithm has complexity $$\\mathcal{O}(N \cdot M^3)\$$. Also, it likely doesn't find the optimal solution unless you use a very small step size, and you have to choose minimum and maximum correctly. There are other algorithms that will both give you a guaranteed optimal solution without having to specify a range and a step size, and which have complexity $$\\mathcal{O}(N)\$$. This is possible since the problem of linear least squares can be solved analytically, and can be implemented efficiently on a computer.

If you don't want to go for the exact mathematical approach, and perhaps want to perform regressions of more complex functions that don't have an exact solution, then you might want to look into algorithms that approach the problem in a smarter way than just brute-forcing all possible values of the parameters. A typical solution to this kind of problem is to implement a gradient descent algorithm.

# Choosing the right error function

The way you calculate the error (or rather the residual to be more precise) is quite strange. Normally one would try to minimize the sum of the squares of the difference between a data point and the predicted value. However, your error function is using the linear difference, and weighting it such that the first datapoint has the least weight, and the last datapoint weighs as much as all the other datapoints combined:

float error = Module(Predict(a, b, c, input[0]) - output[0]);
for (int i = 1; i < input.size(); i++)
error = (Module(Predict(a, b, c, input[i]) - output[i]) + error)/2;


Unless this is really the intended function, you should probably replace this with a function that calculates the least squares residual.

• @Jenia The code you've posted here does not implement gradient descent. Nov 21, 2020 at 15:12