spectral_norm_alioth.js

gistfile1.js

function A(i,j) {
  return 1/((i+j)*(i+j+1)/2+i+1);
}

function Au(u,v) {
  for (var i=0; i<u.length; ++i) {
    var t = 0;
    for (var j=0; j<u.length; ++j)
      t += A(i,j) * u[j];
    v[i] = t;
  }
}

function Atu(u,v) {
  for (var i=0; i<u.length; ++i) {
    var t = 0;
    for (var j=0; j<u.length; ++j)
      t += A(j,i) * u[j];
    v[i] = t;
  }
}

function AtAu(u,v,w) {
  Au(u,w);
  Atu(w,v);
}

function spectralnorm(n) {
  var i, u=[], v=[], w=[], vv=0, vBv=0;
  for (i=0; i<n; ++i) {
    u[i] = 1; v[i] = w[i] = 0; 
  }
  for (i=0; i<10; ++i) {
    AtAu(u,v,w);
    AtAu(v,u,w);
  }
  for (i=0; i<n; ++i) {
    vBv += u[i]*v[i];
    vv  += v[i]*v[i];
  }
  return Math.sqrt(vBv/vv);
}

console.log(spectralnorm(5500).toFixed(9));

spectral_norm.cljs

(ns spectral-norm.core)

(defn a [i j]
  (/ 1 (+ (/ (* (+ i j) (+ i j 1)) 2) i 1)))

(defn mul-av [n v av]
  (dotimes [i n]
    (aset av i 0)
    (dotimes [j n]
      (aset av i (+ (aget av i) (* (a i j) (aget v j)))))))

(defn mul-atv [n v atv]
  (dotimes [i n]
    (aset atv i 0)
    (dotimes [j n]
      (aset atv i (+ (aget atv i) (* (a j i) (aget v j)))))))

(defn mul-atav [n v atav]
  (let [u (array)]
    (mul-av n v u)
    (mul-atv n u atav)))

(defn approximate [n]
  (let [u (array)
        v (array)]
    (dotimes [i n] (aset u i 1) (aset v i 0))
    (dotimes [i 10]
      (mul-atav n u v)
      (mul-atav n v u))
    (loop [i 0 vbv 0 vv 0]
      (if (< i n)
        (recur (inc i)
               (+ vbv (* (aget u i) (aget v i)))
               (+ vv (* (aget v i) (aget v i))))
        (Math/sqrt (/ vbv vv))))))

(defn _print [& args]
  (.log js/console
        (->> args (interpose " ") (apply str))))

(set! *print-fn* _print)

(time
 (.log js/console (approximate 5500)))