Open in app

Sign in

Write

Sign in

Stochastic Gradient Descent in C to factorize recommendation matrices

Kibicho Murage
7 min readSep 27, 2023
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

double learningRate = 0.01;
int maxIterations = 100;
double epsilon = 0.0001;

#define numberOfRows 4
#define numberOfColumns 5
#define middleMultiplier 2

void PrintMatrix(int rowCount, int columnCount, double matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   printf("%.3f ", matrix[i][j]);
  }
  printf("\n");
 }
 printf("\n");
}

void InitializeRandomMatrix(int rowCount, int columnCount, double matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   matrix[i][j] = (double)rand() / RAND_MAX;
  }
 }
}

void MultiplyMatrices(int rowA, int columnA, double A[rowA][columnA], int rowB, int columnB, double B[rowB][columnB], double result[rowA][columnB])
{
 assert(columnA == rowB);
 for(int i = 0; i < rowA; i++)
 {
  for(int j = 0; j < columnB; j++)
  {
   result[i][j] = 0;
   for(int k = 0; k < rowB; k++)
   {
    result[i][j] += A[i][k] * B[k][j];
   }
  }
 }
}

void CalculateError(int rowCount, int columnCount, double A[rowCount][columnCount], double B[rowCount][columnCount], double *error)
{
 *error = 0.0f;
 for(int i = 0; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   *error += pow(A[i][j] - B[i][j] , 2);
  }
 }
}

void FactorizeMatrix(int leftRow, int leftColumn, double left[leftRow][leftColumn],int rightRow, int rightColumn, double right[rightRow][rightColumn],double matrix[leftRow][rightColumn],double matrixResult[leftRow][rightColumn])
{
 assert(leftColumn == rightRow);
 double prediction = 0.0f;
 double error = 0.0f;
 double totalError = 0.0f;
 for(int iteration = 0 ; iteration < maxIterations; iteration++)
 {
  for(int i = 0; i < leftRow; i++)
  {
   for(int j = 0; j < rightColumn; j++)
   {
    //Do non-zero entries
    if(matrix[i][j] != 0)
    {
     prediction = 0.0;
     for(int k = 0; k < leftColumn; k++)
     {
      prediction += left[i][k] * right[k][j];
     }
     error = matrix[i][j] - prediction;
     for(int k = 0; k < leftColumn; k++)
     {
      left[i][k]  += learningRate * (2 * error * right[k][j]);
      right[k][j] += learningRate * (2 * error * left[i][k] );
     }
    }
   }
  }
  MultiplyMatrices(leftRow, leftColumn, left, rightRow, rightColumn, right, matrixResult);
  CalculateError(leftRow, rightColumn, matrix,matrixResult, &totalError);
  if(totalError < epsilon){break;}
 }
}

int main()
{
 double error = 0.0f;
 double matrix[numberOfRows][numberOfColumns] = 
 {
  {5, 0, 1, 0, 4},
  {0, 4, 0, 0, 2},
  {1, 0, 0, 3, 0},
  {0, 2, 3, 0, 0} 
 };
 
 double leftFactor[numberOfRows][middleMultiplier] = {{0.0f}};InitializeRandomMatrix(numberOfRows,middleMultiplier,leftFactor);
 double rightFactor[middleMultiplier][numberOfColumns] = {{0.0f}};InitializeRandomMatrix(middleMultiplier,numberOfColumns,rightFactor);
 double matrixResult[numberOfRows][numberOfColumns] = {{0.0f}};
 
 FactorizeMatrix(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrix,matrixResult);
 MultiplyMatrices(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrixResult);
 PrintMatrix(numberOfRows, numberOfColumns, matrix);
 PrintMatrix(numberOfRows, middleMultiplier, leftFactor);
 PrintMatrix(middleMultiplier, numberOfColumns, rightFactor);
 PrintMatrix(numberOfRows, numberOfColumns, matrixResult);
 //CalculateError(numberOfRows, numberOfColumns, matrix,matrixResult, &error);
 printf("Error: (%.3f)\n", error);
 
 return 0;
}

SGD with GMP — Use the next one. Better memory allocation. Lower leearningRate and maxIterations if seg fault occurs.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <gmp.h>

double learningRate = 0.002;
int maxIterations = 1000;
double epsilon = 0.0001;
int quantization= 0;
mpf_t temporary0;
mpf_t error;
mpf_t totalError;
mpf_t prediction;
mpf_t learningRateHolder;

#define numberOfRows 4
#define numberOfColumns 5
#define middleMultiplier 60

void PrintMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   gmp_printf("%.3Ff ", matrix[i][j]);
  }
  printf("\n");
 }
 printf("\n");
}

void InitializeRandomMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_init(matrix[i][j]);
   mpf_set_d(matrix[i][j], (double)rand() / RAND_MAX);
  }
 }
}

void InitializeZeroMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_init(matrix[i][j]);
   mpf_set_ui(matrix[i][j], 0);
  }
 }
}

void DestroyMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_clear(matrix[i][j]);
  }
 }
}

void MultiplyMatrices(int rowA, int columnA, mpf_t A[rowA][columnA], int rowB, int columnB, mpf_t B[rowB][columnB], mpf_t result[rowA][columnB])
{
 assert(columnA == rowB);
 for(int i = 0; i < rowA; i++)
 {
  for(int j = 0; j < columnB; j++)
  {
   mpf_set_ui(result[i][j], 0);
   for(int k = 0; k < rowB; k++)
   {
    mpf_mul(temporary0, A[i][k], B[k][j]);
    mpf_add(result[i][j], result[i][j], temporary0); 
   }
  }
 }
}


void CalculateError(int rowCount, int columnCount, mpf_t A[rowCount][columnCount], mpf_t B[rowCount][columnCount])
{
 mpf_set_ui(error, 0);
 for(int i = 0; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_sub(temporary0, A[i][j], B[i][j]);
   mpf_mul(temporary0,temporary0,temporary0);
   mpf_add(error,error,temporary0);
  }
 }
}

void QuantizeMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 if(quantization > 0){for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_mul_ui(matrix[i][j],matrix[i][j],quantization);
   mpf_ceil(matrix[i][j], matrix[i][j]);
   mpf_div_ui(matrix[i][j],matrix[i][j],quantization);
  }
 }
 }
}

void FactorizeMatrix(int leftRow, int leftColumn, mpf_t left[leftRow][leftColumn],int rightRow, int rightColumn, mpf_t right[rightRow][rightColumn], mpf_t matrix[leftRow][rightColumn], mpf_t matrixResult[leftRow][rightColumn])
{
 assert(leftColumn == rightRow);
 for(int iteration = 0 ; iteration < maxIterations; iteration++)
 {
  for(int i = 0; i < leftRow; i++)
  {
   for(int j = 0; j < rightColumn; j++)
   {
    //Do non-zero entries
    if(mpf_cmp_ui(matrix[i][j], 0) != 0)
    {
     mpf_set_ui(prediction, 0);
     for(int k = 0; k < leftColumn; k++)
     {
      mpf_mul(temporary0,left[i][k],right[k][j]);
      mpf_add(prediction,prediction,temporary0);
     }
     mpf_sub(error, matrix[i][j], prediction);
     for(int k = 0; k < leftColumn; k++)
     {
      /*Left*/
      mpf_mul(temporary0,error,right[k][j]);
      mpf_mul(temporary0,temporary0, learningRateHolder);
      mpf_add(left[i][k],left[i][k],temporary0);
      /*Right*/
      mpf_mul(temporary0,error,left[i][k]);
      mpf_mul(temporary0,temporary0, learningRateHolder);
      mpf_add(right[k][j],right[k][j],temporary0);
     }
    }
   }
  }
  //MultiplyMatrices(leftRow, leftColumn, left, rightRow, rightColumn, right, matrixResult);
  //CalculateError(leftRow, rightColumn, matrix,matrixResult);
  //if(totalError < epsilon){break;}
 }
 QuantizeMatrix(numberOfRows, middleMultiplier, left);
 QuantizeMatrix(middleMultiplier, numberOfColumns, right);
}

void SetMatrix(int rowCount, int columnCount, mpf_t destination[rowCount][columnCount], double source[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_set_d(destination[i][j], source[i][j]);
  }
 }
}


int main()
{
 quantization = 1000;
 double matrixHolder[numberOfRows][numberOfColumns] = 
 {
  { 12, 15, 17, 18,36},
  { 10, 11, 24, 25, 29},
  { 2,   6,  9, 11, 14},
  { 4,   5,  6, 10, 11} 
 };
 
 
 mpf_init(temporary0);mpf_init(error);mpf_init(totalError);mpf_init(prediction);mpf_init(learningRateHolder);
 mpf_set_d(learningRateHolder, learningRate);
 
 mpf_t matrix[numberOfRows][numberOfColumns];InitializeZeroMatrix(numberOfRows, numberOfColumns, matrix);
 mpf_t leftFactor[numberOfRows][middleMultiplier] = {{0.0f}};InitializeRandomMatrix(numberOfRows,middleMultiplier,leftFactor);
 mpf_t rightFactor[middleMultiplier][numberOfColumns] = {{0.0f}};InitializeRandomMatrix(middleMultiplier,numberOfColumns,rightFactor);
 mpf_t matrixResult[numberOfRows][numberOfColumns] = {{0.0f}};InitializeZeroMatrix(numberOfRows, numberOfColumns, matrixResult);
 
 
 
 SetMatrix(numberOfRows, numberOfColumns, matrix, matrixHolder);
 FactorizeMatrix(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrix,matrixResult);
 MultiplyMatrices(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrixResult);
 PrintMatrix(numberOfRows, numberOfColumns, matrix);
 PrintMatrix(numberOfRows, middleMultiplier, leftFactor);
 PrintMatrix(middleMultiplier, numberOfColumns, rightFactor);
 PrintMatrix(numberOfRows, numberOfColumns, matrixResult);
 CalculateError(numberOfRows, numberOfColumns, matrix,matrixResult);
 gmp_printf("Error: (%.3Ff)\n", error);
 
 DestroyMatrix(numberOfRows, numberOfColumns, matrix);
 DestroyMatrix(numberOfRows, middleMultiplier, leftFactor);
 DestroyMatrix(middleMultiplier, numberOfColumns, rightFactor);
 mpf_clear(temporary0);mpf_clear(error);mpf_clear(totalError);mpf_clear(prediction);mpf_clear(learningRateHolder);
 return 0;
}
//Fixed Memory allocations
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <gmp.h>
//https://stackoverflow.com/questions/74848782/problem-of-memory-allocation-after-operations-on-2d-gmp-integer-matrix
double learningRate = 0.0002;
int maxIterations = 400;
double epsilon = 0.0001;
int quantization= 0;
mpf_t temporary0;
mpf_t error;
mpf_t totalError;
mpf_t prediction;
mpf_t learningRateHolder;

#define numberOfRows 18
#define numberOfColumns 40
#define middleMultiplier 1

void PrintMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   gmp_printf("%6.1Ff  ", matrix[i][j]);
  }
  printf("\n");
 }
 printf("\n");
}

void SetRandomMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_set_d(matrix[i][j], (double)rand() / RAND_MAX);
  }
 }
}

void SetZeroMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_set_ui(matrix[i][j], 0);
  }
 }
}

void InitializeMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_init(matrix[i][j]);
  }
 }
}



void DestroyMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_clear(matrix[i][j]);
  }
 }
 free(matrix);
}

void MultiplyMatrices(int rowA, int columnA, mpf_t A[rowA][columnA], int rowB, int columnB, mpf_t B[rowB][columnB], mpf_t result[rowA][columnB])
{
 assert(columnA == rowB);
 for(int i = 0; i < rowA; i++)
 {
  for(int j = 0; j < columnB; j++)
  {
   mpf_set_ui(result[i][j], 0);
   for(int k = 0; k < rowB; k++)
   {
    mpf_mul(temporary0, A[i][k], B[k][j]);
    mpf_add(result[i][j], result[i][j], temporary0); 
   }
  }
 }
}


void CalculateError(int rowCount, int columnCount, mpf_t A[rowCount][columnCount], mpf_t B[rowCount][columnCount])
{
 mpf_set_ui(error, 0);
 for(int i = 0; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_sub(temporary0, A[i][j], B[i][j]);
   mpf_mul(temporary0,temporary0,temporary0);
   mpf_add(error,error,temporary0);
  }
 }
}

void QuantizeMatrix(int rowCount, int columnCount, mpf_t matrix[rowCount][columnCount])
{
 if(quantization > 0){for(int i = 0 ; i < rowCount; i++)
 {
  for(int j = 0; j < columnCount; j++)
  {
   mpf_mul_ui(matrix[i][j],matrix[i][j],quantization);
   mpf_ceil(matrix[i][j], matrix[i][j]);
   mpf_div_ui(matrix[i][j],matrix[i][j],quantization);
  }
 }
 }
}

void FactorizeMatrix(int leftRow, int leftColumn, mpf_t left[leftRow][leftColumn],int rightRow, int rightColumn, mpf_t right[rightRow][rightColumn], mpf_t matrix[leftRow][rightColumn], mpf_t matrixResult[leftRow][rightColumn])
{
 assert(leftColumn == rightRow);
 for(int iteration = 0 ; iteration < maxIterations; iteration++)
 {
  for(int i = 0; i < leftRow; i++)
  {
   for(int j = 0; j < rightColumn; j++)
   {
    //Do non-zero entries
    if(mpf_cmp_ui(matrix[i][j], 0) != 0)
    {
     mpf_set_ui(prediction, 0);
     for(int k = 0; k < leftColumn; k++)
     {
      mpf_mul(temporary0,left[i][k],right[k][j]);
      mpf_add(prediction,prediction,temporary0);
     }
     mpf_sub(error, matrix[i][j], prediction);
     for(int k = 0; k < leftColumn; k++)
     {
      /*Left*/
      mpf_mul(temporary0,error,right[k][j]);
      mpf_mul(temporary0,temporary0, learningRateHolder);
      mpf_add(left[i][k],left[i][k],temporary0);
      /*Right*/
      mpf_mul(temporary0,error,left[i][k]);
      mpf_mul(temporary0,temporary0, learningRateHolder);
      mpf_add(right[k][j],right[k][j],temporary0);
     }
    }
   }
  }
  //MultiplyMatrices(leftRow, leftColumn, left, rightRow, rightColumn, right, matrixResult);
  //CalculateError(leftRow, rightColumn, matrix,matrixResult);
  //if(totalError < epsilon){break;}
 }
 QuantizeMatrix(numberOfRows, middleMultiplier, left);
 QuantizeMatrix(middleMultiplier, numberOfColumns, right);
}

void SetMatrix(int rowCount, int columnCount, mpf_t destination[rowCount][columnCount], double source[rowCount][columnCount])
{
 double sum = 0;
 for(int i = 0 ; i < rowCount; i++)
 {
  sum = 0;
  for(int j = 0; j < columnCount; j++)
  {
   sum += source[i][j]/10;
   mpf_set_d(destination[i][j], sum);
   //mpf_set_d(destination[i][j], source[i][j]);
  }
 }
}

void EncodeFile(char *fileName)
{
 int c = 0; int count = 0; int i = 0; int j = 0;
 FILE *fp = fopen(fileName, "rb");
 assert(fp != NULL);
 
 mpf_init(temporary0);mpf_init(error);mpf_init(totalError);mpf_init(prediction);mpf_init(learningRateHolder);
 mpf_set_d(learningRateHolder, learningRate);

 mpf_t (*matrix)[numberOfRows] = malloc( sizeof(mpf_t[numberOfRows][numberOfColumns]));InitializeMatrix(numberOfRows, numberOfColumns, matrix);

 mpf_t (*leftFactor)[numberOfRows] = malloc( sizeof(mpf_t[numberOfRows][middleMultiplier]));InitializeMatrix(numberOfRows,middleMultiplier,leftFactor);
 mpf_t (*rightFactor)[middleMultiplier] = malloc( sizeof(mpf_t[middleMultiplier][numberOfColumns]));InitializeMatrix(middleMultiplier,numberOfColumns,rightFactor);
 mpf_t (*matrixResult)[numberOfRows] = malloc( sizeof(mpf_t[numberOfRows][numberOfColumns]));InitializeMatrix(numberOfRows, numberOfColumns, matrixResult);
 
 double matrixHolder[numberOfRows][numberOfColumns];
 quantization = 10;
 while((c = fgetc(fp)) != EOF)
 {
  //printf("(%d %d)%d| ",i,j,c);
  matrixHolder[i][j] = (double) c;
  j+=1;
  if(j % numberOfColumns == 0)
  {
   j = 0;
   i += 1;
   if(i % numberOfRows == 0)
   {
    printf("%d:\n",count);
    SetMatrix(numberOfRows, numberOfColumns, matrix, matrixHolder);
    SetRandomMatrix(numberOfRows, middleMultiplier, leftFactor);
    SetRandomMatrix(middleMultiplier, numberOfColumns, rightFactor);
    
    FactorizeMatrix(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrix,matrixResult);
    MultiplyMatrices(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrixResult);
    PrintMatrix(numberOfRows, middleMultiplier, leftFactor);
    PrintMatrix(middleMultiplier, numberOfColumns, rightFactor);
    PrintMatrix(numberOfRows, numberOfColumns, matrixResult);
    PrintMatrix(numberOfRows, numberOfColumns, matrix);
    CalculateError(numberOfRows, numberOfColumns, matrix,matrixResult);
    
    if(count == 100){break;}
    i = 0;j = 0;count+=1;
   }
  }
 }
 fclose(fp);
 DestroyMatrix(numberOfRows, numberOfColumns, matrix);
 DestroyMatrix(numberOfRows, middleMultiplier, leftFactor);
 DestroyMatrix(middleMultiplier, numberOfColumns, rightFactor);
 DestroyMatrix(numberOfRows, numberOfColumns, matrixResult);
 mpf_clear(temporary0);mpf_clear(error);mpf_clear(totalError);mpf_clear(prediction);mpf_clear(learningRateHolder);
}

int main()
{
 char *fileName = "cat2.mp4";
 EncodeFile(fileName);
 return 0;
}
/*
int main()
{
 quantization = 1000;
 char *fileName = "cat2.mp4";
 
 double matrixHolder[numberOfRows][numberOfColumns] = 
 {
  { 12, 15, 17, 18,36},
  { 10, 11, 24, 25, 29},
  { 2,   6,  9, 11, 14},
  { 4,   5,  6, 10, 11} 
 };
 
 
 mpf_init(temporary0);mpf_init(error);mpf_init(totalError);mpf_init(prediction);mpf_init(learningRateHolder);
 mpf_set_d(learningRateHolder, learningRate);
 
 mpf_t matrix[numberOfRows][numberOfColumns];InitializeZeroMatrix(numberOfRows, numberOfColumns, matrix);
 mpf_t leftFactor[numberOfRows][middleMultiplier] = {{0.0f}};InitializeRandomMatrix(numberOfRows,middleMultiplier,leftFactor);
 mpf_t rightFactor[middleMultiplier][numberOfColumns] = {{0.0f}};InitializeRandomMatrix(middleMultiplier,numberOfColumns,rightFactor);
 mpf_t matrixResult[numberOfRows][numberOfColumns] = {{0.0f}};InitializeZeroMatrix(numberOfRows, numberOfColumns, matrixResult);
 
 
 
 SetMatrix(numberOfRows, numberOfColumns, matrix, matrixHolder);
 FactorizeMatrix(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrix,matrixResult);
 MultiplyMatrices(numberOfRows, middleMultiplier, leftFactor, middleMultiplier, numberOfColumns, rightFactor, matrixResult);
 PrintMatrix(numberOfRows, numberOfColumns, matrix);
 PrintMatrix(numberOfRows, middleMultiplier, leftFactor);
 PrintMatrix(middleMultiplier, numberOfColumns, rightFactor);
 PrintMatrix(numberOfRows, numberOfColumns, matrixResult);
 CalculateError(numberOfRows, numberOfColumns, matrix,matrixResult);
 gmp_printf("Error: (%.3Ff)\n", error);
 
 DestroyMatrix(numberOfRows, numberOfColumns, matrix);
 DestroyMatrix(numberOfRows, middleMultiplier, leftFactor);
 DestroyMatrix(middleMultiplier, numberOfColumns, rightFactor);
 mpf_clear(temporary0);mpf_clear(error);mpf_clear(totalError);mpf_clear(prediction);mpf_clear(learningRateHolder);
 return 0;
}*/

--

--

Kibicho Murage

Written by Kibicho Murage

0 Followers

AI Researcher at Fileforma. Twitter : murage_kibicho

Help

Status

About

Careers

Blog

Privacy

Terms

Text to speech

Teams