feat: add ml/strided/dsgd-trainer-squared-epsilon-insensitive

lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/examples/index.js

@@ -0,0 +1,61 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2024 The Stdlib Authors.
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*    http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+'use strict';
+
+var uniform = require( '@stdlib/random/array/uniform' );
+var Float64Array = require( '@stdlib/array/float64' );
+var dsgdTrainerSqEpsIns = require( './../lib/base.js' );
+
+var M = 100000; // number of samples
+var N = 20;     // number of features
+
+var opts = {
+	'dtype': 'float64'
+};
+
+// Generate a random `MxN` design matrix and a random target vector:
+var x = uniform( M*N, -10.0, 10.0, opts );
+var y = uniform( M, -50.0, 50.0, opts );
+
+// Allocate the weight vector and the per-feature L1 shrinkage workspace:
+var w = new Float64Array( N );
+var workspace = new Float64Array( N );
+
+// Configure the training hyperparameters:
+var penalty = 'elasticnet';
+var fitIntercept = true;
+var l1Ratio = 0.15;
+var maxIter = 20;
+var learningRate = 'invscaling';
+var eta0 = 0.01;
+var powerT = 0.25;
+var epsilon = 0.0;
+var lambda = 1.0e-5;
+var intercept = 0.0;
+
+// Train a linear model via plain SGD with an elastic-net penalty:
+var out = dsgdTrainerSqEpsIns( penalty, learningRate, fitIntercept, M, N, l1Ratio, maxIter, eta0, powerT, epsilon, lambda, intercept, y, 1, 0, w, 1, 0, x, N, 1, 0, workspace, 1, 0 ); // eslint-disable-line max-len
+
+var j;
+console.log( 'Estimated intercept: %d', out.intercept );
+console.log( '' );
+console.log( 'feature | estimated weight' );
+for ( j = 0; j < N; j++ ) {
+	console.log( '%d\t| %d', j, out.weights[ j ] );
+}

lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/lib/base.js

@@ -0,0 +1,256 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2024 The Stdlib Authors.
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*    http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+'use strict';
+
+// MODULES //
+
+var logger = require( 'debug' );
+var sqEpsInsGradient = require( '@stdlib/ml/base/loss/float64/squared-epsilon-insensitive-gradient' );
+var ddot = require( '@stdlib/blas/base/ddot' ).ndarray;
+var daxpy = require( '@stdlib/blas/base/daxpy' ).ndarray;
+var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
+var pow = require( '@stdlib/math/base/special/pow' );
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+
+
+// VARIABLES //
+
+var debug = logger( 'ml:dsgd-trainer-squared-epsilon-insensitive' );
+
+var MIN_SCALE = 1.0e-9;
+var MAX_DLOSS = 1e12;
+
+// NOTE: both of the above constants were picked from sklearn API
+
+// FUNCTIONS //
+
+/**
+* Applies the cumulative L1 penalty to the weight vector using scikit-learn's truncated-gradient scheme.
+*
+* @private
+* @param {NonNegativeInteger} N - number of features (length of `w`)
+* @param {number} u - cumulative L1 penalty accumulated so far
+* @param {number} scaleFactor - current scaling factor applied to the stored weights
+* @param {Float64Array} w - weight vector (updated in-place)
+* @param {integer} strideW - stride length for `w`
+* @param {NonNegativeInteger} offsetW - starting index for `w`
+* @param {Float64Array} workspace - workspace array tracking the total L1 shrinkage applied to each feature (updated in-place)
+* @param {integer} strideWS - stride length for `workspace`
+* @param {NonNegativeInteger} offsetWS - starting index for `workspace`
+* @returns {void}
+*/
+function l1Penalty( N, u, scaleFactor, w, strideW, offsetW, workspace, strideWS, offsetWS ) {
+	var wsIdx;
+	var idx;
+	var z;
+	var j;
+
+	for ( j = 0; j < N; j++ ) {
+		idx = offsetW + ( j * strideW );
+		wsIdx = offsetWS + ( j * strideWS );
+		z = w[ idx ];
+		if ( scaleFactor * z > 0.0 ) {
+			w[ idx ] = max( 0.0, z - ( ( u + workspace[ wsIdx ] ) / scaleFactor ) );
+		} else if ( scaleFactor * z < 0.0 ) {
+			w[ idx ] = min( 0.0, z + ( ( u - workspace[ wsIdx ] ) / scaleFactor ) );
+		}
+		workspace[ wsIdx ] += scaleFactor * ( w[ idx ] - z );
+	}
+}
+
+
+// MAIN //
+
+/**
+* Trains a linear model with the squared epsilon-insensitive loss via stochastic gradient descent.
+*
+* @private
+* @param {string} penalty - penalty type: `'l1'`, `'l2'`, or `'elasticnet'`
+* @param {string} learningRate - schedule: `'constant'` (eta = eta0), `'invscaling'` (eta = eta0/t^powerT), `'basic'` (eta = 10/(10+t)), or `'pegasos'` (eta = 1/(lambda*t))
+* @param {boolean} fitIntercept - boolean indicating whether to fit an intercept term
+* @param {NonNegativeInteger} M - number of samples
+* @param {NonNegativeInteger} N - number of features
+* @param {number} l1Ratio - elastic-net mixing parameter on the interval `[0,1]` (`0` => L2, `1` => L1); overridden for pure L1/L2 penalties
+* @param {PositiveInteger} maxIter - number of epochs
+* @param {number} eta0 - initial/base learning rate (used by `'constant'` and `'invscaling'`)
+* @param {number} powerT - exponent for the inverse-scaling schedule (only used by `'invscaling'`)
+* @param {number} epsilon - width of the insensitive region of the loss
+* @param {number} lambda - regularization strength
+* @param {number} intercept - initial intercept
+* @param {Float64Array} y - target vector
+* @param {integer} strideY - stride length for `y`
+* @param {NonNegativeInteger} offsetY - starting index for `y`
+* @param {Float64Array} w - weight vector
+* @param {integer} strideW - stride length for `w`
+* @param {NonNegativeInteger} offsetW - starting index for `w`
+* @param {Float64Array} x - `M` by `N` input matrix
+* @param {integer} strideX1 - stride of the first dimension of `x`
+* @param {integer} strideX2 - stride of the second dimension of `x`
+* @param {NonNegativeInteger} offsetX - starting index for `x`
+* @param {Float64Array} workspace - workspace array
+* @param {integer} strideWS - stride length for `workspace`
+* @param {NonNegativeInteger} offsetWS - starting index for `workspace`
+* @returns {Object} results object with `intercept` and `weights`
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* // Two features, four samples; target y = 2*x0 - 1*x1:
+* var x = new Float64Array( [ 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 1.0 ] );
+* var y = new Float64Array( [ 2.0, -1.0, 1.0, 3.0 ] );
+* var w = new Float64Array( 2 );
+* var workspace = new Float64Array( 2 );
+*
+* var out = dsgdTrainerSqEpsIns(
+*     'l2', 'invscaling', true, // penalty, learningRate, fitIntercept
+*     4, 2, 0.0,                // M, N, l1Ratio
+*     1000,                     // maxIter
+*     0.02, 0.5,                // eta0, powerT
+*     0.0, 1.0e-4,              // epsilon, lambda
+*     0.0,                      // initial intercept
+*     y, 1, 0,                  // y, strideY, offsetY
+*     w, 1, 0,                  // w, strideW, offsetW
+*     x, 2, 1, 0,               // x, strideX1, strideX2, offsetX
+*     workspace, 1, 0           // workspace, strideWS, offsetWS
+* );
+*/
+function dsgdTrainerSqEpsIns( penalty, learningRate, fitIntercept, M, N, l1Ratio, maxIter, eta0, powerT, epsilon, lambda, intercept, y, strideY, offsetY, w, strideW, offsetW, x, strideX1, strideX2, offsetX, workspace, strideWS, offsetWS ) { // eslint-disable-line max-params, max-len
+	var scaleFactor;
+	var update;
+	var factor;
+	var epoch;
+	var dloss;
+	var eta;
+	var ox;
+	var oy;
+	var t;
+	var p;
+	var i;
+	var u;
+
+	debug( 'Starting SGD trainer with squared epsilon insensitive loss...' );
+	debug( 'M = %d, N = %d, penalty = %s, lr scheduler = %s, loss function = %s, eta0 = %s, ', M, N, penalty, learningRate, 'squared-epsilon-insensitive-loss', eta0 );
+
+	/*
+	* w - (k*N) (feature vector)
+	* x - M*N (input vector)
+	* y - k (output vector)
+	*/
+
+	/*
+	* DOUBTS?
+	* - should I include the intercept with the weight vector itself?
+	*    - if yes, then should N be including intercept?
+	*    - if no, then how do I show it to the user? should I have a `intercept` field in the results struct
+	*
+	*    - Should we implement a standalone weight matrix just like how `ml/incr/sgd-regression` as well as
+	*        the sklearn API did? They referred sofia-ml as far as I know.
+	*/
+
+	// Do we set this here or expect higher level user to pass it?
+	if ( penalty === 'l2' ) {
+		l1Ratio = 0.0;
+	} else if ( penalty === 'l1' ) {
+		l1Ratio = 1.0;
+	}
+
+	eta = eta0;
+	t = 1;
+	scaleFactor = 1.0;
+	u = 0.0;
+	for ( epoch = 1; epoch <= maxIter; epoch++ ) {
+		ox = offsetX;
+		oy = offsetY;
+		for ( i = 0; i < M; i++ ) {
+			p = ( scaleFactor*ddot( N, w, strideW, offsetW, x, strideX2, ox ) ) + intercept; // eslint-disable-line max-len
+
+			if ( learningRate === 'invscaling' ) {
+				eta = eta0 / pow( t, powerT );
+			} else if ( learningRate === 'basic' ) {
+				eta = 10.0 / ( 10.0 + t );
+			} else if ( learningRate === 'pegasos' ) {
+				eta = 1.0 / ( lambda * t );
+			}
+			// NOTE: No need to write branch for `constant` as we set eta=eta0 before starting iterations
+
+			dloss = sqEpsInsGradient( 1.0, epsilon, y[ oy ], p );
+			if ( dloss < -MAX_DLOSS ) {
+				dloss = -MAX_DLOSS;
+			} else if ( dloss > MAX_DLOSS ) {
+				dloss = MAX_DLOSS;
+			}
+			update = -eta * dloss;
+
+			if ( penalty === 'l2' || penalty === 'elasticnet' ) {
+				factor = 1.0 - ( ( 1.0 - l1Ratio ) * eta * lambda );
+				scaleFactor *= max( 0.0, factor );
+			}
+
+			if ( scaleFactor < MIN_SCALE ) {
+				dscal( N, scaleFactor, w, strideW, offsetW );
+				scaleFactor = 1.0;
+			}
+
+			// Gradient step: w_eff += -eta*dloss*x  <=>  w_stored += (-eta*dloss/scaleFactor)*x
+			// w_eff = w_stored*scaleFactor
+			daxpy( N, update/scaleFactor, x, strideX2, ox, w, strideW, offsetW ); // eslint-disable-line max-len
+
+			if ( fitIntercept ) {
+				intercept += update;
+			}
+
+			if ( penalty === 'l1' || penalty === 'elasticnet' ) {
+				u += ( l1Ratio * eta * lambda );
+				l1Penalty( N, u, scaleFactor, w, strideW, offsetW, workspace, strideWS, offsetWS ); // eslint-disable-line max-len
+			}
+
+			t += 1;
+			ox += strideX1;
+			oy += strideY;
+		}
+		debug( 'Epoch %d done...', epoch );
+	}
+
+	// Fold any residual scale back into the weights so `w` holds effective values:
+	// sklearn does w.reset_scale()
+	if ( scaleFactor !== 1 ) {
+		dscal( N, scaleFactor, w, strideW, offsetW );
+		scaleFactor = 1.0;
+	}
+
+	debug( 'Finished SGD trainer with squared epsilon insensitive loss.' );
+
+	return {
+		'intercept': intercept,
+		'weights': w
+	};
+
+	/*
+	* returning this doesn't look solid, I have to brainstorm here.
+	* we will have to either return `intercept`, or keep `intercept` inside the weight matrix,
+	* as we cant pass a number by reference as fn argument, so the `intercept` user passes won't be updated, instead
+	* a copy of that `intercept` is being updated here.
+	*/
+}
+
+
+// EXPORTS //
+
+module.exports = dsgdTrainerSqEpsIns;