Public Member Functions
void	initialize_Z (parameters &weights)

void	compute_AZx ()

void	update_eigenvalues ()

void	compute_delta ()

void	update_K ()

void	update_A ()

void	update_b ()

void	update_D ()

void	check ()

	~OjaNewton ()

Public Attributes
vw *	all

std::shared_ptr< rand_state >	_random_state

int	m

int	epoch_size

float	alpha

int	cnt

int	t

float *	ev

float *	b

float *	D

float **	A

float **	K

float *	zv

float *	vv

float *	tmp

example **	buffer

float *	weight_buffer

struct update_data	data

float	learning_rate_cnt

bool	normalize

bool	random_init

Detailed Description

Definition at line 32 of file OjaNewton.cc.

Constructor & Destructor Documentation

◆ ~OjaNewton()

OjaNewton::~OjaNewton ( )

inline

Definition at line 345 of file OjaNewton.cc.

   {
     free(ev);
     free(b);
     free(D);
     free(buffer);
     free(weight_buffer);
     free(zv);
     free(vv);
     free(tmp);
     if (A)
     {
       for (int i = 1; i <= m; i++)
       {
         free(A[i]);
         free(K[i]);
       }
     }
 
     free(A);
     free(K);
 
     free(data.Zx);
     free(data.AZx);
     free(data.delta);
   }

Member Function Documentation

◆ check()

void OjaNewton::check ( )

inline

Definition at line 264 of file OjaNewton.cc.

References vw::num_bits, parameters::strided_index(), and vw::weights.

Referenced by learn().

   {
     double max_norm = 0;
     for (int i = 1; i <= m; i++)
       for (int j = i; j <= m; j++) max_norm = fmax(max_norm, fabs(K[i][j]));
     // printf("|K| = %f\n", max_norm);
     if (max_norm < 1e7)
       return;
 
     // implicit -> explicit representation
     // printf("begin conversion: t = %d, norm(K) = %f\n", t, max_norm);
 
     // first step: K <- AKA'
 
     // K <- AK
     for (int j = 1; j <= m; j++)
     {
       memset(tmp, 0, sizeof(double) * (m + 1));
 
       for (int i = 1; i <= m; i++)
       {
         for (int h = 1; h <= m; h++)
         {
           tmp[i] += A[i][h] * K[h][j];
         }
       }
 
       for (int i = 1; i <= m; i++) K[i][j] = tmp[i];
     }
     // K <- KA'
     for (int i = 1; i <= m; i++)
     {
       memset(tmp, 0, sizeof(double) * (m + 1));
 
       for (int j = 1; j <= m; j++)
         for (int h = 1; h <= m; h++) tmp[j] += K[i][h] * A[j][h];
 
       for (int j = 1; j <= m; j++)
       {
         K[i][j] = tmp[j];
       }
     }
 
     // second step: w[0] <- w[0] + (DZ)'b, b <- 0.
 
     uint32_t length = 1 << all->num_bits;
     for (uint32_t i = 0; i < length; i++)
     {
       weight& w = all->weights.strided_index(i);
       for (int j = 1; j <= m; j++) w += (&w)[j] * b[j] * D[j];
     }
 
     memset(b, 0, sizeof(double) * (m + 1));
 
     // third step: Z <- ADZ, A, D <- Identity
 
     // double norm = 0;
     for (uint32_t i = 0; i < length; ++i)
     {
       memset(tmp, 0, sizeof(float) * (m + 1));
       weight& w = all->weights.strided_index(i);
       for (int j = 1; j <= m; j++)
       {
         for (int h = 1; h <= m; ++h) tmp[j] += A[j][h] * D[h] * (&w)[h];
       }
       for (int j = 1; j <= m; ++j)
       {
         // norm = std::max(norm, fabs(tmp[j]));
         (&w)[j] = tmp[j];
       }
     }
     // printf("|Z| = %f\n", norm);
 
     for (int i = 1; i <= m; i++)
     {
       memset(A[i], 0, sizeof(double) * (m + 1));
       D[i] = 1;
       A[i][i] = 1;
     }
   }

◆ compute_AZx()

void OjaNewton::compute_AZx ( )

inline

Definition at line 117 of file OjaNewton.cc.

Referenced by learn().

   {
     for (int i = 1; i <= m; i++)
     {
       data.AZx[i] = 0;
       for (int j = 1; j <= i; j++)
       {
         data.AZx[i] += A[i][j] * data.Zx[j];
       }
     }
   }

◆ compute_delta()

void OjaNewton::compute_delta ( )

inline

Definition at line 147 of file OjaNewton.cc.

References f.

Referenced by learn().

   {
     data.bdelta = 0;
     for (int i = 1; i <= m; i++)
     {
       float gamma = fmin(learning_rate_cnt / t, 1.f);
 
       // if different learning rates are used
       /*data.delta[i] = gamma * data.AZx[i] * data.sketch_cnt;
       for (int j = 1; j < i; j++) {
           data.delta[i] -= A[i][j] * data.delta[j];
       }
       data.delta[i] /= A[i][i];*/
 
       // if a same learning rate is used
       data.delta[i] = gamma * data.Zx[i] * data.sketch_cnt;
 
       data.bdelta += data.delta[i] * b[i];
     }
   }

◆ initialize_Z()

void OjaNewton::initialize_Z ( parameters & weights )

inline

Definition at line 60 of file OjaNewton.cc.

References f, NORM2, vw::num_bits, and parameters::strided_index().

Referenced by save_load().

   {
     uint32_t length = 1 << all->num_bits;
     if (normalize)  // initialize normalization part
     {
       for (uint32_t i = 0; i < length; i++) (&(weights.strided_index(i)))[NORM2] = 0.1f;
     }
     if (!random_init)
     {
       // simple initialization
       for (int i = 1; i <= m; i++) (&(weights.strided_index(i)))[i] = 1.f;
     }
     else
     {
       // more complicated initialization: orthgonal basis of a random matrix
 
       const double PI2 = 2.f * 3.1415927f;
 
       for (uint32_t i = 0; i < length; i++)
       {
         weight& w = weights.strided_index(i);
         float r1, r2;
         for (int j = 1; j <= m; j++)
         {
           // box-muller tranform: https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
           // redraw until r1 should be strictly positive
           do
           {
             r1 = _random_state->get_and_update_random();
             r2 = _random_state->get_and_update_random();
           } while (r1 == 0.f);
 
           (&w)[j] = std::sqrt(-2.f * log(r1)) * (float)cos(PI2 * r2);
         }
       }
     }
 
     // Gram-Schmidt
     for (int j = 1; j <= m; j++)
     {
       for (int k = 1; k <= j - 1; k++)
       {
         double tmp = 0;
 
         for (uint32_t i = 0; i < length; i++)
           tmp += ((double)(&(weights.strided_index(i)))[j]) * (&(weights.strided_index(i)))[k];
         for (uint32_t i = 0; i < length; i++)
           (&(weights.strided_index(i)))[j] -= (float)tmp * (&(weights.strided_index(i)))[k];
       }
       double norm = 0;
       for (uint32_t i = 0; i < length; i++)
         norm += ((double)(&(weights.strided_index(i)))[j]) * (&(weights.strided_index(i)))[j];
       norm = std::sqrt(norm);
       for (uint32_t i = 0; i < length; i++) (&(weights.strided_index(i)))[j] /= (float)norm;
     }
   }

◆ update_A()

void OjaNewton::update_A ( )

inline

Definition at line 182 of file OjaNewton.cc.

Referenced by learn().

   {
     for (int i = 1; i <= m; i++)
     {
       for (int j = 1; j < i; j++)
       {
         zv[j] = 0;
         for (int k = 1; k <= i; k++)
         {
           zv[j] += A[i][k] * K[k][j];
         }
       }
 
       for (int j = 1; j < i; j++)
       {
         vv[j] = 0;
         for (int k = 1; k <= j; k++)
         {
           vv[j] += A[j][k] * zv[k];
         }
       }
 
       for (int j = 1; j < i; j++)
       {
         for (int k = j; k < i; k++)
         {
           A[i][j] -= vv[k] * A[k][j];
         }
       }
 
       float norm = 0;
       for (int j = 1; j <= i; j++)
       {
         float temp = 0;
         for (int k = 1; k <= i; k++)
         {
           temp += K[j][k] * A[i][k];
         }
         norm += A[i][j] * temp;
       }
       norm = sqrtf(norm);
 
       for (int j = 1; j <= i; j++)
       {
         A[i][j] /= norm;
       }
     }
   }

◆ update_b()

void OjaNewton::update_b ( )

inline

Definition at line 231 of file OjaNewton.cc.

Referenced by learn().

   {
     for (int j = 1; j <= m; j++)
     {
       float tmp = 0;
       for (int i = j; i <= m; i++)
       {
         tmp += ev[i] * data.AZx[i] * A[i][j] / (alpha * (alpha + ev[i]));
       }
       b[j] += tmp * data.g;
     }
   }

◆ update_D()

void OjaNewton::update_D ( )

inline

Definition at line 244 of file OjaNewton.cc.

   {
     for (int j = 1; j <= m; j++)
     {
       float scale = fabs(A[j][j]);
       for (int i = j + 1; i <= m; i++) scale = fmin(fabs(A[i][j]), scale);
       if (scale < 1e-10)
         continue;
       for (int i = 1; i <= m; i++)
       {
         A[i][j] /= scale;
         K[j][i] *= scale;
         K[i][j] *= scale;
       }
       b[j] /= scale;
       D[j] *= scale;
       // printf("D[%d] = %f\n", j, D[j]);
     }
   }

◆ update_eigenvalues()

void OjaNewton::update_eigenvalues ( )

inline

Definition at line 129 of file OjaNewton.cc.

References f.

Referenced by learn().

   {
     for (int i = 1; i <= m; i++)
     {
       float gamma = fmin(learning_rate_cnt / t, 1.f);
       float tmp = data.AZx[i] * data.sketch_cnt;
 
       if (t == 1)
       {
         ev[i] = gamma * tmp * tmp;
       }
       else
       {
         ev[i] = (1 - gamma) * t * ev[i] / (t - 1) + gamma * t * tmp * tmp;
       }
     }
   }

◆ update_K()

void OjaNewton::update_K ( )

inline

Definition at line 168 of file OjaNewton.cc.

Referenced by learn().

   {
     float tmp = data.norm2_x * data.sketch_cnt * data.sketch_cnt;
     for (int i = 1; i <= m; i++)
     {
       for (int j = 1; j <= m; j++)
       {
         K[i][j] += data.delta[i] * data.Zx[j] * data.sketch_cnt;
         K[i][j] += data.delta[j] * data.Zx[i] * data.sketch_cnt;
         K[i][j] += data.delta[i] * data.delta[j] * tmp;
       }
     }
   }