@@ -246,3 +246,147 @@ export class TSNE {
246246 }
247247}
248248
249+ export class MDS {
250+ nComponents : number ;
251+ metric : boolean ;
252+ nInit : number ;
253+ maxIter : number ;
254+ eps : number ;
255+
256+ embedding_ : Float64Array [ ] | null = null ;
257+ stress_ : number | null = null ;
258+
259+ constructor (
260+ options : {
261+ nComponents ?: number ;
262+ metric ?: boolean ;
263+ nInit ?: number ;
264+ maxIter ?: number ;
265+ eps ?: number ;
266+ } = { } ,
267+ ) {
268+ this . nComponents = options . nComponents ?? 2 ;
269+ this . metric = options . metric ?? true ;
270+ this . nInit = options . nInit ?? 4 ;
271+ this . maxIter = options . maxIter ?? 300 ;
272+ this . eps = options . eps ?? 1e-3 ;
273+ }
274+
275+ fitTransform ( X : Float64Array [ ] ) : Float64Array [ ] {
276+ const n = X . length ;
277+ // Compute distance matrix
278+ const D = new Float64Array ( n * n ) ;
279+ for ( let i = 0 ; i < n ; i ++ ) {
280+ for ( let j = i + 1 ; j < n ; j ++ ) {
281+ let d = 0 ;
282+ const xi = X [ i ] ?? new Float64Array ( 0 ) ;
283+ const xj = X [ j ] ?? new Float64Array ( 0 ) ;
284+ for ( let k = 0 ; k < xi . length ; k ++ )
285+ d += ( ( xi [ k ] ?? 0 ) - ( xj [ k ] ?? 0 ) ) ** 2 ;
286+ d = Math . sqrt ( d ) ;
287+ D [ i * n + j ] = d ;
288+ D [ j * n + i ] = d ;
289+ }
290+ }
291+
292+ // Classical MDS via double centering
293+ const d = this . nComponents ;
294+ // B = -0.5 * H * D^2 * H where H = I - (1/n) * 11^T
295+ const D2 = new Float64Array ( n * n ) ;
296+ for ( let i = 0 ; i < n * n ; i ++ ) D2 [ i ] = ( D [ i ] ?? 0 ) ** 2 ;
297+
298+ const rowMean = new Float64Array ( n ) ;
299+ const colMean = new Float64Array ( n ) ;
300+ let totalMean = 0 ;
301+ for ( let i = 0 ; i < n ; i ++ ) {
302+ for ( let j = 0 ; j < n ; j ++ ) {
303+ rowMean [ i ] = ( rowMean [ i ] ?? 0 ) + ( D2 [ i * n + j ] ?? 0 ) ;
304+ colMean [ j ] = ( colMean [ j ] ?? 0 ) + ( D2 [ i * n + j ] ?? 0 ) ;
305+ totalMean += D2 [ i * n + j ] ?? 0 ;
306+ }
307+ }
308+ for ( let i = 0 ; i < n ; i ++ ) {
309+ rowMean [ i ] = ( rowMean [ i ] ?? 0 ) / n ;
310+ colMean [ i ] = ( colMean [ i ] ?? 0 ) / n ;
311+ }
312+ totalMean /= n * n ;
313+
314+ const B = new Float64Array ( n * n ) ;
315+ for ( let i = 0 ; i < n ; i ++ ) {
316+ for ( let j = 0 ; j < n ; j ++ ) {
317+ B [ i * n + j ] =
318+ - 0.5 *
319+ ( ( D2 [ i * n + j ] ?? 0 ) -
320+ ( rowMean [ i ] ?? 0 ) -
321+ ( colMean [ j ] ?? 0 ) +
322+ totalMean ) ;
323+ }
324+ }
325+
326+ // Power iteration to get top-d eigenvectors of B
327+ const vecs : Float64Array [ ] = [ ] ;
328+ const vals : number [ ] = [ ] ;
329+ const Bcopy = new Float64Array ( B ) ;
330+ for ( let comp = 0 ; comp < d ; comp ++ ) {
331+ const v = new Float64Array ( n ) ;
332+ for ( let i = 0 ; i < n ; i ++ ) v [ i ] = Math . random ( ) - 0.5 ;
333+ for ( let iter = 0 ; iter < 100 ; iter ++ ) {
334+ const w = new Float64Array ( n ) ;
335+ for ( let i = 0 ; i < n ; i ++ ) {
336+ for ( let j = 0 ; j < n ; j ++ )
337+ w [ i ] ! += ( Bcopy [ i * n + j ] ?? 0 ) * ( v [ j ] ?? 0 ) ;
338+ }
339+ let norm = 0 ;
340+ for ( let i = 0 ; i < n ; i ++ ) norm += ( w [ i ] ?? 0 ) ** 2 ;
341+ norm = Math . sqrt ( norm ) || 1 ;
342+ for ( let i = 0 ; i < n ; i ++ ) v [ i ] = ( w [ i ] ?? 0 ) / norm ;
343+ if ( iter === 99 ) {
344+ let lam = 0 ;
345+ for ( let i = 0 ; i < n ; i ++ ) lam += ( w [ i ] ?? 0 ) * ( v [ i ] ?? 0 ) ;
346+ vals . push ( lam ) ;
347+ }
348+ }
349+ vecs . push ( v ) ;
350+ // Deflate
351+ const lam = vals [ comp ] ?? 0 ;
352+ for ( let i = 0 ; i < n ; i ++ ) {
353+ for ( let j = 0 ; j < n ; j ++ ) {
354+ Bcopy [ i * n + j ] ! -= lam * ( v [ i ] ?? 0 ) * ( v [ j ] ?? 0 ) ;
355+ }
356+ }
357+ }
358+
359+ // Embedding: X_new[i][k] = sqrt(lambda_k) * v_k[i]
360+ const Y : Float64Array [ ] = Array . from (
361+ { length : n } ,
362+ ( ) => new Float64Array ( d ) ,
363+ ) ;
364+ for ( let k = 0 ; k < d ; k ++ ) {
365+ const scale = Math . sqrt ( Math . max ( vals [ k ] ?? 0 , 0 ) ) ;
366+ for ( let i = 0 ; i < n ; i ++ ) {
367+ ( Y [ i ] as Float64Array ) [ k ] = scale * ( ( vecs [ k ] as Float64Array ) [ i ] ?? 0 ) ;
368+ }
369+ }
370+
371+ this . embedding_ = Y ;
372+ // Compute stress
373+ let stress = 0 ;
374+ for ( let i = 0 ; i < n ; i ++ ) {
375+ for ( let j = i + 1 ; j < n ; j ++ ) {
376+ let distY = 0 ;
377+ const yi = Y [ i ] as Float64Array ;
378+ const yj = Y [ j ] as Float64Array ;
379+ for ( let k = 0 ; k < d ; k ++ ) distY += ( ( yi [ k ] ?? 0 ) - ( yj [ k ] ?? 0 ) ) ** 2 ;
380+ distY = Math . sqrt ( distY ) ;
381+ stress += ( distY - ( D [ i * n + j ] ?? 0 ) ) ** 2 ;
382+ }
383+ }
384+ this . stress_ = stress ;
385+ return Y ;
386+ }
387+
388+ fit ( X : Float64Array [ ] ) : this {
389+ this . fitTransform ( X ) ;
390+ return this ;
391+ }
392+ }
0 commit comments