diff --git a/lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/examples/index.js b/lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/examples/index.js new file mode 100644 index 000000000000..94fd74d8e05a --- /dev/null +++ b/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 ] ); +} diff --git a/lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/lib/base.js b/lib/node_modules/@stdlib/ml/strided/dsgd-trainer-squared-epsilon-insensitive/lib/base.js new file mode 100644 index 000000000000..20796fbfafa5 --- /dev/null +++ b/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;