2011年6月26日,星期日

更快的LDA:第二部分

我能够提高其他效率 Vowpal兔子的LDA实施 通过 向量化。基本思想是找到重复的昂贵操作的循环并通过 SIMD 同时执行多个昂贵的操作。特别是我正在利用 上证所 指令系统。

向量化范例

对于我的简单功能 近似对数和指数,转换为SSE是一项最直接的机械任务(如此机械,实际上,我想知道为什么编译器无法为我完成此任务?)。这是一个示例:首先,原始的快速以2为底的对数,
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

inline float
fastlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  float y = vx.i;
  y *= 1.0f / (1 << 23);

  return
    y - 124.22544637f - 1.498030302f * mx.f - 1.72587999f / (0.3520887068f + mx.f);
}
接下来是SSE版本
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

typedef float v4sf __attribute__ ((__vector_size__ (16)));
typedef int v4si __attribute__ ((__vector_size__ (16)));
#define v4sf_dup_literal(x) ((v4sf) { x, x, x, x })
#define v4si_dup_literal(x) ((v4si) { x, x, x, x })
#define v4sf_from_v4si __builtin_ia32_cvtdq2ps

inline v4sf
vfastlog2 (v4sf x)
{
  union { v4sf f; v4si i; } vx = { x };
  union { v4si i; v4sf f; } mx = {
    (vx.i & v4si_dup_literal (0x007FFFFF)) | v4si_dup_literal ((0x7e << 23))
  };
  v4sf y = v4sf_from_v4si (vx.i);
  y *= v4sf_dup_literal ((1.0f / (1 << 23)));
  
  const v4sf c_124_22544637 = v4sf_dup_literal (124.22544637f);
  const v4sf c_1_498030302 = v4sf_dup_literal (1.498030302f);
  const v4sf c_1_725877999 = v4sf_dup_literal (1.72587999f);
  const v4sf c_0_3520087068 = v4sf_dup_literal (0.3520887068f);
  
  return 
    y - c_124_22544637 - c_1_498030302 * mx.f - c_1_725877999 / (c_0_3520087068 + mx.f);
}
vfastlog2 运行速度与 fastlog2 (实际上有时会稍快一些,具体取决于 fastlog2 并使用x87运算进行编译),并计算几乎相同的内容(有时差异很小,比近似误差小得多)。但是,它可以一次执行4次计算。

有一个陷阱:它总是一次执行4次计算,因此,如果数组的长度不是16个字节的倍数,则最后一次计算将从内存的其他位读取(可能导致分段错误)。此外,输入必须对齐16字节,否则会出现分段错误。有人告诉我,在超高性能应用程序中,开发人员总是将阵列排列为16字节对齐,并且是16字节长的倍数。但是,没有简单的方法来安排此操作,因此我不得不进行一些防御性的编码。我使用的模式是
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

typedef float v4sf_aligned __attribute__ ((__vector_size__ (16))) __attribute__ ((aligned (16)));

float f (float);
v4sf vectorized_f (v4sf);

inline void
conservative_map (float* x,
                  size_t n)
{
  unsigned int i;

  // handle non-aligned prefix
  for (i = 0; i < n && (((uintptr_t) (x + i)) % 16) > 0; ++i)
    { 
      x[i] = f (x[i]);
    }

  // the good stuff
  for (; i + 4 < n; i += 4)
    { 
      * (v4sf_aligned *) (x + i) = vectorized_f (* (v4sf_aligned*) (x + i));
    }

  // handle remainder
  for (; i < n; ++i)
    { 
      x[i] = f (x[i]);
    }
}
对于大约100个主题的地图,这在最坏的情况下会导致25个操作之外的6个附加操作,因此平均可能会有10%的罚款。

结果

我在与我的基准相同的基准上进行了一些其他计时测试 以前的帖子。 \ [
\ begin {array} {c | c}
\ mbox {代码版本}&\ mbox {总挂钟} \\ \ hline
\ mbox {原始}&\ mbox {79秒} \\
\ mbox {近似值,无SSE}&\ mbox {49秒} \\
\ mbox {近似值,含SSE}&\ mbox {39秒}
\ end {array}
\]
近似加速比矢量化加速大,这很好,因为即使对于不支持SSE的环境,近似加速始终可用。

2011年6月24日,星期五

更快的LDA

现在,我已经掌握了 字节重新解释的邪恶力量 在整数和IEEE浮点数之间,我一直在寻找加速的东西。最近发布的 快速实施LDA 由Smola等等提醒我看Vowpal Wabbit的 LDA实施。的 官方网站 该文件充满了看起来昂贵的先验功能,可以使用近似版本来加快速度,希望不会破坏学习算法。

更快的Log Gamma和Digamma

除了许多我已经有了快速近似的指数和对数调用之外,还有一些对 对数伽玛地伽玛。幸运的是, 斯特林近似 是对数,因此我可以利用我的快速对数近似值 以前的博客文章.

我发现对于Log Gamma和Digamma而言,两次应用移位公式可以在精度和速度之间取得最佳折衷。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

inline float
fastlgamma (float x)
{
  float logterm = fastlog (x * (1.0f + x) * (2.0f + x));
  float xp3 = 3.0f + x;

  return
    -2.081061466f - x + 0.0833333f / xp3 - logterm + (2.5f + x) * fastlog (xp3);
}

inline float
fastdigamma (float x)
{
  float twopx = 2.0f + x;
  float logterm = fastlog (twopx);

  return - (1.0f + 2.0f * x) / (x * (1.0f + x))
         - (13.0f + 6.0f * x) / (12.0f * twopx * twopx)
         + logterm;
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$ x \ in [1/100,10] $的平均相对准确度(使用 lgammaf 来自libc作为Log Gamma的黄金标准,并且 boost :: math :: digamma 作为Digamma的黄金标准)。
fastlgamma relative accuracy = 0.00045967
fastdigamma relative accuracy = 0.000420604
fastlgamma million calls per second = 60.259
lgammaf million calls per second = 21.4703
boost::math::lgamma million calls per second = 1.8951
fastdigamma million calls per second = 66.0639
boost::math::digamma million calls per second = 18.9672
好吧,Log Gamma的速度大约是3倍(但是为什么boost :: math :: lgamma这么慢?)却略好于Digamma的3倍。

LDA时序

所以我将这些近似值放入 官方网站 并使用 wiki1K 来自的文件 测试/训练集/ 目录;但是我将文件本身连接了100次以减慢速度。
原始运行
% (cd test; time ../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatch 128 train-sets/wiki1Kx100.dat)
your learning rate is too high, setting it to 1
using no cache
Reading from train-sets/wiki1Kx100.dat
num sources = 1
Num weight bits = 13
learning rate = 1
initial_t = 1
power_t = 0.5
learning_rate set to 1
average    since       example  example    current  current  current
loss       last        counter   weight      label  predict features
10.296875  10.296875         3      3.0    unknown   0.0000       38
10.437156  10.577436         6      6.0    unknown   0.0000       14
10.347227  10.239314        11     11.0    unknown   0.0000       32
10.498633  10.650038        22     22.0    unknown   0.0000        2
10.495566  10.492500        44     44.0    unknown   0.0000      166
10.469184  10.442189        87     87.0    unknown   0.0000       29
10.068007  9.666831        174    174.0    unknown   0.0000       17
9.477440   8.886873        348    348.0    unknown   0.0000        2
9.020482   8.563524        696    696.0    unknown   0.0000      143
8.482232   7.943982       1392   1392.0    unknown   0.0000       31
8.106654   7.731076       2784   2784.0    unknown   0.0000       21
7.850269   7.593883       5568   5568.0    unknown   0.0000       25
7.707427   7.564560      11135  11135.0    unknown   0.0000       39
7.627417   7.547399      22269  22269.0    unknown   0.0000       61
7.583820   7.540222      44537  44537.0    unknown   0.0000        5
7.558824   7.533827      89073  89073.0    unknown   0.0000      457

finished run
number of examples = 0
weighted example sum = 0
weighted label sum = 0
average loss = -nan
best constant = -nan
total feature number = 0
../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatc  69.06s user 13.60s system 101% cpu 1:21.59 total
近似运行
% (cd test; time ../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatch 128 train-sets/wiki1Kx100.dat)               your learning rate is too high, setting it to 1
using no cache
Reading from train-sets/wiki1Kx100.dat
num sources = 1
Num weight bits = 13
learning rate = 1
initial_t = 1
power_t = 0.5
learning_rate set to 1
average    since       example  example    current  current  current
loss       last        counter   weight      label  predict features
10.297077  10.297077         3      3.0    unknown   0.0000       38
10.437259  10.577440         6      6.0    unknown   0.0000       14
10.347348  10.239455        11     11.0    unknown   0.0000       32
10.498796  10.650243        22     22.0    unknown   0.0000        2
10.495748  10.492700        44     44.0    unknown   0.0000      166
10.469374  10.442388        87     87.0    unknown   0.0000       29
10.068179  9.666985        174    174.0    unknown   0.0000       17
9.477574   8.886968        348    348.0    unknown   0.0000        2
9.020435   8.563297        696    696.0    unknown   0.0000      143
8.482178   7.943921       1392   1392.0    unknown   0.0000       31
8.106636   7.731093       2784   2784.0    unknown   0.0000       21
7.849980   7.593324       5568   5568.0    unknown   0.0000       25
7.707124   7.564242      11135  11135.0    unknown   0.0000       39
7.627207   7.547283      22269  22269.0    unknown   0.0000       61
7.583952   7.540696      44537  44537.0    unknown   0.0000        5
7.559260   7.534566      89073  89073.0    unknown   0.0000      457

finished run
number of examples = 0
weighted example sum = 0
weighted label sum = 0
average loss = -nan
best constant = -nan
total feature number = 0
../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatc  41.97s user 7.69s system 102% cpu 48.622 total
因此,从81秒降低到49秒,速度提高了大约40%。但这有效吗?好吧,我让马特(Matt)担任最终裁判(我给他寄了一封修改过的 官方网站 (用于测试),但初步迹象表明,使用昂贵功能的近似版本不会破坏学习过程。

下一步:矢量化

肯定会有更多的速度,因为 官方网站 加载像
for (size_t i = 0; i<global.lda; i++)
    {
      sum += gamma[i];
      gamma[i] = mydigamma(gamma[i]);
    }
只是为了尖叫 向量化。敬请关注!

2011年6月23日,星期四

减少AUC到逐点回归

查尔斯·埃尔坎洛杉矶机器学习聚会昨天。在演讲之前,我们正在讨论AUC损失的排名问题。众所周知 归结为成对分类 导致AUC后悔界限为$ 2 r $,其中$ r $是对归类分类问题的遗憾。查尔斯说,在实践中,他通过简化为逐点逻辑回归获得了很好的结果,即估计``相关概率''并使用归纳排名。他推测后勤损失实际上是在某种程度上限制了AUC损失,这是因为它给数据点提供了有效的重要性权重(高表示``令人惊讶'',低表示``不令人惊讶'');有趣的是,它比平方损失要好得多;铰链丢失根本不起作用。关于铰链损耗的观察是有道理的,因为已知减少逐点分类并不可靠。但是平方或逻辑损失呢?

因此,在考虑减少量时,存在一个一致性问题:对引起的子问题零后悔会导致对原始问题的零后悔吗?由于平方和逻辑损失是 正确的评分规则,一致性在直觉上是可能的。但是,对于嘈杂的问题,还存在效率问题。换句话说,对引发问题的后悔如何限制对原始问题的后悔?

物流损失是棘手的,但是平方损失我取得了一些进展。假设从$ S = X \ times Y $中抽取示例对,其中$ Y = \ {0,1 \} $根据分布$ D_ {S \ times} = D_x \ times D_ {y | x} \ times D_x \ times D_ {y | x} $。如果您熟悉排名问题,那么您可能正在调整,因为在实践中很少有明确的逐点相关性二进制概念。但是,由于我正在谈论归结为点回归的减少,因此我将假设它可以取得进展。

我现在要修复$ x_1 $和$ x_2 $,最后我会期望$ D_ {x | 0} \ times D_ {x | 1} $。给定分类器$ h:S \ to \ mathbb {R} $,条件AUC损失由\ [
\ mathrm {AUC}(h | x_1,x_2)= E _ {(y_1,y_2)\ sim D _ {(y_1,y_2)| (x_1,x_2)}} \ left [\ left(1_ {y_1> y_2} 1_{h (x_1) < h (x_2)} + 1_{y_1 < y_2} 1_{h (x_1) >h(x_2)} \ right)\ right]。
\]同时平方损失由\ [
L_2(h | x_1,x_2)= E _ {(y_1,y_2)\ sim D _ {(y_1,y_2)| (x_1,x_2)}} \ left [\ sum_k y_k(1-h(x_k))^ 2 +(1-y_k)h(x_k)^ 2 \ right]。
\]
不失一般性地假定$ h(x_1)<h(x_2)$。然后有条件的AUC损失变为\ [
\ mathrm {AUC}(h_1,h_2,p_1,p_2)= p_1(1-p_2),
\]其中$ h_k = h(x_k)$和$ p_k = p(y_k = 1 | x_k)$。最佳AUC损失由\ [
\ mathrm {AUC}(h_1,h_2,p_1,p_2)= \ begin {cases}
p_1(1-p_2)和\ mbox {if} p_1 \ leq p_2,\\
(1-p_1)p_2& \mbox{if } p_1 > p_2.
\ end {cases}
\],因此AUC后悔由\ [
r_ \ mathrm {AUC} = \ max \ left(p_1-p_2,0 \ right)。
\]同时对分类器的平方损失遗憾是\ [
\ begin {aligned}
r_2(h; x_1,x_2)&= \sum_k r_2 (h_k;x_k)= \ sum_k(p_k-h_k)^ 2。
\ end {aligned}
\]作为对手,我的目标是针对给定的平方损失后悔最大化AUC后悔。要解决的程序是\ [\ max_ {p_1,p_2,h_1,h_2}(p_1-p_2)\]受\ [
\ begin {aligned}
p_2-p_1&\ leq 0,\\
h_1-h_2&\ leq 0,\\
\ sum_k(p_k-h_k)^ 2&= r_2 (h; x_1, x_2).
\ end {aligned}
\]通过 KKT条件 (感谢Mathematica!)产生解\ [
\ begin {aligned}
h_1&= h_2,\\
h_2&= \ frac {p_1 + p_2} {2},\\
p_1-p_2&= \sqrt{2 r_2 (h; x_1, x_2)}.
\ end {aligned}
\]这就是说,对手的最佳选择是将$ h_1 $和$ h_2 $放置在$ p_1 $和$ p_2 $的均值之上和之下的$ \ epsilon $;这让人联想到 中位数选民定理。在这种情况下,AUC的遗憾是\ [
r _ {\ mathrm {AUC}}(h; x_1,x_2)\ leq \ sqrt {2 r_2(h | x_1,x_2)}。
\]现在我可以对两边的$ D_ {x | 0} \ times D_ {x | 1} $做出期望,\ [
\ begin {aligned}
r _ {\ mathrm {AUC}}(h)&\leq E_{(x_1, x_2) \sim D_{x|0} \times D_{x|1}} \left[ \sqrt{2 r_2 (h;x_1,x_2)} \ right] \\
&\leq \sqrt{ 2 E_{(x_1, x_2) \sim D_{x|0} \times D_{x|1}} \left[ r_2(h; x_1,x_2)\right] } \\
&= \sqrt{2 \left( E_{x \sim D_{x|0}} \left[ r (h;x)\ right] + E_ {x \ sim D_ {x | 1}} \ left [r(h; x)\ right)\ right)} \\
&= 2 \sqrt{E_{x \sim \tilde D_x} \left[ r (h; x) \right]} \\
&= 2 \ sqrt {\ tilde r_2(h)}
\ end {aligned}
\]这里$ D_ {x | y} $是给定标签$ y $的$ x $的条件分布,$ \ tilde D_x = \ frac {D_ {x | 0} + D_ {x | 1}} {2 } $是$ x $的类平衡分布,而$ \ tilde r_2(h)$是基础分类器相对于$ \ tilde D_x $的平方误差后悔。那么这是什么意思?
  • 如果针对类平衡的数据集(即相同数目的正负数)计算平方损失后悔,则此分析仅保证减少平方损失的一致性。这可以通过例如重要性加权或拒绝采样来实现。不平衡的数据集可能会出现问题,例如,当使用几乎完全是负样本的数据集和一个正样本的数据集时,基础数据集的平均每样本均方差损失后悔与AUC后悔规模之间的界限为数据集(因为在一个积极的例子上的重要性越来越大,但被忽略了)。
  • 后悔界限对潜在遗憾具有平方根依赖性,这比成对归类到分类时对潜在遗憾的线性依赖性更差。

您可能会想到,后勤损失不适合该技术。但是,我开始欣赏查尔斯发表评论的直觉。想象一下一个在线学习场景,一个物流学习者在不考虑他们的班级标签的情况下就被不断地输入示例。如果数据集非常不平衡,例如学习者大多是负面的,他们将很快收敛于一个高度负面的常数。这将有效降低负面示例的学习率,而正面示例将保持较高的学习率。这模仿了对数据加权以进行类平衡的重要性操作。如何将这种直觉转化为正式声明?

2011年6月20日,星期一

快速近似对数,指数,幂和逆根

受启发 伊恩·斯蒂芬森(Ian Stephenson),我发现了一个更快更准确的$ \ log_2 $近似值,这反过来又提高了我对$ p $根的反向估计。所产生的技巧通常对于计算能力$ x ^ p $运作良好,这也是计算$ L_p $范数的一部分。还有一个快速的近似对数和指数,它们在广义上相关,因为这些运算在实践中被广泛地用于 表示概率.

更新资料:这里的代码示例已过期;小问题已修复,您应该从 Google代码存储库。 (更新资料2:Google代码已失效,但是 版本存在于github

快速近似对数

快速的\ log_2 $近似值意味着对自然对数或编译时已知的任何其他基数的快速近似值,因为差是一个乘法。因此,考虑一个浮点​​数$ x $和值$ I_x $,它们将$ x $的字节表示形式解释为整数\ [
\ begin {aligned}
x&=(1 + m_x)2 ^ {e_x},\\
\ log_2(x)&= e_x + \ log_2(1 + m_x),\\
I_x&=(e_x + B)L + m_x L,\\
I_x / L-B&= e_x + m_x,\\
\ log_2(x)&= I_x / L-B + \ log_2(1 + m_x)-m_x。
\ end {aligned}
\]因此,需要一种从$ x $中提取尾数并近似为[0,1)$中的$ z \ g_2(1 + z)-z $的尾数的方法。我采用了有理近似值,这似乎是准确性和速度的最佳结合点,
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  float y = vx.i;
  y *= 1.0 / (1 << 23);

  return
    y - 124.22544637f - 1.498030302f * mx.f - 1.72587999f / (0.3520887068f + mx.f);
}

inline float
fastlog (float x)
{
  return 0.69314718f * fastlog2 (x);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {-2 + \ log(8)} {\ log(4)} $
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  float y = vx.i;
  y *= 1.0 / (1 << 23);
  return y - 126.94269504f;
}

inline float
fasterlog (float x)
{
  return 0.69314718f * fasterlog2 (x);
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$ x \ in [1/100,10] $中的平均相对精度。
fastlog2 relative accuracy = 2.09352e-05
fastlog relative accuracy = 2.09348e-05
fasterlog2 relative accuracy = 0.0130367
fasterlog relative accuracy = 0.0130367
fastlog2 million calls per second = 160.141
fastlog million calls per second = 143.552
fasterlog2 million calls per second = 218.345
fasterlog million calls per second = 210.435
log2f million calls per second = 40.8511

快速近似指数

同样,2是最方便的基础。这里$ 2 ^ x $是我要近似的值,$ I_ {2 ^ x} $是与将$ 2 ^ x $的字节模式解释为整数相关的值。 \ [
\ begin {aligned}
2 ^ x&= 2 ^ {x-\ lfloor x \ rfloor} 2 ^ {\ lfloor x \ rfloor} \\
&=(1 + 2 ^ {x-\ lfloor x \ rfloor}-1)2 ^ {\ lfloor x \ rfloor},\\
I_ {2 ^ x}&=(\ lfloor x \ rfloor + B)L +(2 ^ {x-\ lfloor x \ rfloor}-1)L \\
&= L(x + B-1 + 2 ^ {x-\ lfloor x \ rfloor}-x + \ lfloor x \ rfloor)。
\ end {aligned}
\]校正项的格式为$ g(z)= 2 ^ z-z $ for $ z \ in [0,1)$,并且我又一次从有理函数逼近中得到了最大的收益。这导致
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastpow2 (float p)
{
  union { float f; uint32_t i; } vp = { p };
  int sign = (vp.i >> 31);
  int w = p;
  float z = p - w + sign;
  union { uint32_t i; float f; } v = { (1 << 23) * (p + 121.2740838f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z) };
  return v.f;
}

inline float
fastexp (float p)
{
  return fastpow2 (1.442695040f * p);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {2-\ log(2)} {2 \ log(2)} $
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterpow2 (float p)
{
  union { uint32_t i; float f; } v = { (1 << 23) * (p + 126.94269504f) };
  return v.f;
}

inline float
fasterexp (float p)
{
  return fasterpow2 (1.442695040f * p);
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还在[1/20,20] $中测量了$ p的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow2 relative accuracy (positive p) = 1.58868e-05
fastexp relative accuracy (positive p) = 1.60712e-05
fasterpow2 relative accuracy (positive p) = 0.0152579
fasterexp relative accuracy (positive p) = 0.0152574
fastpow2 relative accuracy (inverse root p) = 1.43517e-05
fastexp relative accuracy (inverse root p) = 1.7255e-05
fasterpow2 relative accuracy (inverse root p) = 0.013501
fasterexp relative accuracy (inverse root p) = 0.0111832
fastpow2 million calls per second = 153.561
fastexp million calls per second = 143.311
fasterpow2 million calls per second = 215.006
fasterexp million calls per second = 214.44
expf million calls per second = 4.16527

快速近似功率

快速近似对数和快速近似指数的组合产生快速近似幂。
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastpow (float x,
         float p)
{
  return fastpow2 (p * fastlog2 (x));
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow relative accuracy (positive p) = 0.000165618
fastpow relative accuracy (inverse root p) = 0.00011997
fastpow million calls per second = 58.8533
powf million calls per second = 8.44577

快速近似逆根

通过反根,我的意思是为$ p \ geq 1 $解决$ y = x ^ {-1 / p} $。这与快速近似功率不一样吗?是的,但是在这种情况下,以下方法相当准确且速度更快。初始近似值基于修改后的 反平方根hack,\ [
I_y =-\ frac {1} {p} I_x + \ frac {p + 1} {p}(B-\ sigma)L。
\]使用近似精炼 哈雷法 在$ f(y)= \ log(y)+ \ frac {1} {p} \ log(x)$,\ [
\ begin {aligned}
y_ {n + 1}&= y_n \ frac {1-\ delta} {1 + \ delta} \ approx y_n(1-\ delta)^ 2 \\
\ delta&= \ frac {\ log(2)} {2} \ left(\ log_2(y)+ \ frac {1} {p} \ log_2(x)\ right)
\ end {aligned}
\] $ \ log_2 $是使用上面的快速近似对数计算的。总共变成
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterinvproot (float x,
                float p)
{
  union { float f; uint32_t i; } v = { x };
  unsigned int R = 0x5F375A84U;   // R = (3/2) (B - sigma) L
  unsigned int pR = ((2 * (p + 1)) / (3 * p)) * R;
  unsigned int sub = v.i / p;
  v.i = pR - sub;
  return v.f;
}

inline float
fastinvproot (float x,
              float p)
{
  float y = fasterinvproot (x, p);
  float log2x = fastlog2 (x);
  float log2y = fastlog2 (y);
  float err = 1.0 - 0.34657359f * (log2x / p + log2y);
  return y * err * err;
}
注意: fastinvproot 是改进(更快,更准确)的 fast_inverse_p_one 呈现在我的 以前的博客文章.
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastinvproot relative accuracy (positive p) = 0.000727901
fastinvproot relative accuracy (inverse root p) = 0.00300208
fastinvproot million calls per second = 82.7019
所以这有点令人困惑,但是这里的``正p''情况意味着要取一个反$ p $的根,而``反根p''的意思实际上是试图让这个例程陷入计算$ y = x ^ p $。因此,这些数字表示(出于我不完全理解的原因) fastinvproot 查找逆根比使用一般幂函数更准确。

买者自负

除了精度较低外,这些例程还不能处理浮点数 特殊价值 像标准例程一样。另外,如果您不适当地向这些例程提供负输入,则会发生奇怪的事情。

2011年6月16日,星期四

快速逆根

最近,我想知道如何快速近似一个数字的$ p $根,\ [
y = x ^ {-\ frac {1} {p}},
\]使用时出现 $ L_p $归一化.

当$ p = 2 $时 反平方根hack 适用,因此我在Wikipedia上对$ p \ neq 2 $修改了派生结果,导致\ [
I_y =-\ frac {1} {p} I_x + \ frac {p + 1} {p}(B-\ sigma)L,
\]导致了以下代码,
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p (float x,
                float p)
{
  union { float f; uint32_t i; } v = { x };

  unsigned int R = 0x5F375A84U;   // R = (3/2) (B - sigma) L
  unsigned int pR = ((2 * (p + 1)) / (3 * p)) * R;
  unsigned int sub = v.i / p;

  v.i = pR - sub;

  return v.f;
}
这个例程做的很合理,但是缺乏 牛顿法 反平方根hack具有的抛光步骤选项。由于要避免调用powf,因此我将牛顿方法应用于$ f(y)= \ log(y)+ \ frac {1} {p} \ log(x)$,这样会产生迭代\ [
y_ {n + 1} = y_n \ left(1-\ log(2)\ left(\ frac {1} {p} \ log_2(x)+ \ log_2(y)\ right)\ right)。
\]的想法是将其与$ \ log_2 $的廉价近似值一起使用,以利用浮点表示形式。我尝试的第一个近似是一个有理近似,得出了以下结论。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p_one (float x,
                    float p)
{
  float y = fast_inverse_p (x, p);
  union { float f; uint32_t i; } vx = { x };
  union { float f; uint32_t i; } vy = { y };

  // equivalent to: mx = frexp (x, &ex), ex += 126
  int ex = ((vx.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  int ey = ((vy.i >> 23) & 0xFF);
  union { uint32_t i; float f; } my = { (vy.i & 0x007FFFFF) | (0x7e << 23) };

  float log2x = ex - 123.52964741f - 4.29802897f / (0.73958333f + mx.f);
  float log2y = ey - 123.52964741f - 4.29802897f / (0.73958333f + my.f);

  return y * (1.0f - 0.69314718f * (log2x / p + log2y));
}
那有些精确,有些贵。为了获得更高的精度,我采用了多项式逼近法和两次牛顿迭代法。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p_two (float x,
                    float p)
{
  float y = fast_inverse_p (x, p);
  union { float f; uint32_t i; } vx = { x };
  union { float f; uint32_t i; } vy = { y };

  int ex = ((vx.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  int ey = ((vy.i >> 23) & 0xFF);
  union { uint32_t i; float f; } my = { (vy.i & 0x007FFFFF) | (0x7e << 23) };

  float powmx = mx.f;
  float log2x = ex - 129.78533457081f + 10.07763081376f * powmx;
  powmx *= mx.f; log2x -= 13.90982855008f * powmx;
  powmx *= mx.f; log2x += 12.64852864324f * powmx;
  powmx *= mx.f; log2x -= 6.39532347130f * powmx;
  powmx *= mx.f; log2x += 1.364335673877f * powmx;
  log2x /= p;

  float powmy = my.f;
  float log2y = ey - 129.78533457081f + 10.07763081376f * powmy;
  powmy *= my.f; log2y -= 13.90982855008f * powmy;
  powmy *= my.f; log2y += 12.64852864324f * powmy;
  powmy *= my.f; log2y -= 6.39532347130f * powmy;
  powmy *= my.f; log2y += 1.364335673877f * powmy;

  union { float f; uint32_t i; } vz = { y * (1.0f - 0.69314718f * (log2x + log2y)) };
  int ez = ((vz.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mz = { (vz.i & 0x007FFFFF) | (0x7e << 23) };

  float powmz = mz.f;
  float log2z = ez - 129.78533457081f + 10.07763081376f * powmz;
  powmz *= mz.f; log2z -= 13.90982855008f * powmz;
  powmz *= mz.f; log2z += 12.64852864324f * powmz;
  powmz *= mz.f; log2z -= 6.39532347130f * powmz;
  powmz *= mz.f; log2z += 1.364335673877f * powmz;

  return vz.f * (1.0f - 0.69314718f * (log2x + log2z));
}
这更准确,更昂贵。

这是一些示例值以及一些相对详尽的输入范围内的平均相对误差估计。
x       p       fast0   fast1   fast2   powf
7.000   0.870   0.10889 0.10663 0.10681 0.10681
7.000   2.488   0.45789 0.45704 0.45744 0.45744
7.000   4.106   0.63713 0.62284 0.62256 0.62256
7.000   5.724   0.73334 0.71151 0.71181 0.71181
7.000   7.342   0.78715 0.76660 0.76718 0.76718
7.000   8.960   0.82151 0.80411 0.80480 0.80479
69.000  0.870   0.00749 0.00769 0.00770 0.00770
69.000  2.488   0.18674 0.18229 0.18236 0.18236
69.000  4.106   0.36593 0.35650 0.35659 0.35658
69.000  5.724   0.47131 0.47700 0.47726 0.47726
69.000  7.342   0.56050 0.56189 0.56176 0.56176
69.000  8.960   0.63580 0.62382 0.62341 0.62341
211.000 0.870   0.00217 0.00213 0.00213 0.00213
211.000 2.488   0.11642 0.11628 0.11636 0.11636
211.000 4.106   0.27032 0.27150 0.27161 0.27161
211.000 5.724   0.40273 0.39225 0.39260 0.39260
211.000 7.342   0.47678 0.48220 0.48243 0.48243
211.000 8.960   0.54817 0.55025 0.55031 0.55030
err0 = 0.021138
err1 = 0.000680451
err2 = 7.20003e-06
我做了一个时序测试,用 -O3 -finline功能-fast-math,在运行Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。
fast_inverse_p million calls per second = 206.677
fast_inverse_p_one million calls per second = 77.9861
fast_inverse_p_two million calls per second = 26.531
powf million calls per second = 9.04992
巧合的是,每次精度提高都会使吞吐量损失大约三倍。

2011年6月13日,星期一

更好的标签相似性可视化

好吧,我花了整个周末的大部分时间来改善自己的 标签相似度可视化,没有很好的理由,只是发痒时就必须挠痒。这就是最终看起来最好的东西。
  1. 计算每个哈希标签的词项频率向量。
    1. 对于每个推文X:
      1. 对于X $中的每个主题标签$ h:
        1. 对于X $中的每个非标签令牌$ t \:
          1. $(h,t)$的增量计数。
  2. 计算前1000个最频繁的标签中的每个标签对的成对余弦$ \ cos \ theta(h,h ^ \ prime)$。
  3. 将相似度定义为$ s(h,h ^ \ prime)= \ arccos(\ cos \ theta(h,h ^ \ prime))$。
    1. 与$ 1.0-\ cos \ theta(h,h ^ \ prime)$相比,我得到的负特征值更少。这是在超球面上移动和跨越超弦之间的区别。
  4. 做60维 MDS 在$ s $。
    1. 60个特征值中有两个为负,因此我将其视为0。所以实际上是58维MDS。
    2. 我在这里得到任何负特征值时感到有些惊讶,因为我所有的项向量都占据了超球面的正超象限。显然,我的超直觉需要超调……或者我有一个错误。
  5. 将产生的60维表示输入 吨位.
    1. 我使用困惑10。
我选择t-SNE是因为高维结构似乎是相对孤立的簇。我推断是因为我将邻域定义为$ \ epsilon \ theta $球并运行 弗洛伊德·沃歇尔 在邻域图中,我注意到我必须使用一个非常大的邻域($ \ cos \ theta \ geq 0.8 $),然后才能获得足够大的连接组件以变得有趣。

最终,当我绘制此标签时,我尝试将颜色随机化,以便在所有标签相互叠放时能够看到一些东西。真的png不能公正,您应该得到 pdf版本 并放大。

2011年6月9日,星期四

标签相似性可视化

永远不要低估漂亮图片的力量。一方面, 探索性数据分析 确实有助于建立更好的模型。此外,无论您想要更多的补助金还是更多的薪水,漂亮的图片都是说服人们您正在做有趣事情的好方法。

所以我开始研究Twitter 井号 采用;最终目标是通过在编写推文时自动建议或纠正哈希标签,在阅读推文时解释哈希标签等来改进我们的客户端软件。Twitter的新手用户经常抱怨这没有任何意义,因此这将是一个不错的选择增值。

标签受累于 词汇问题。基本上,您说的是#potato,我说的是#spud。通常,人们使用的变体在语义上几乎是等效的,但有时在词汇上却完全不同。一些例子是
  • #shoutouts和#shoutout:这种情况似乎适合词法分析。
  • #imjustsaying和#imjustsayin:这些都非常接近,也许我们希望基于``阻止''的策略能够奏效。
  • #nowplaying和#np:有很多用长缩写代替长哈希标签的示例,以节省费用,也许可以为此制定定制的词汇策略。
  • #libya和#feb17:完全不是词法,但是这些标签的当前几乎等效的性质可能随时间而不稳定。
代替词法驱动的方法,最好开发一种数据驱动的方法来建议,更正或解释哈希标签:毕竟,不乏推文。

好了,偷偷摸摸,这是漂亮的照片。我在推文样本中提取了1000个最受欢迎的哈希标签,并通过将所有带有该哈希标签的推文相加来计算每个哈希标签的词项频率向量。接下来,我计算 余弦相似度 哈希标记之间,阈值设为0.985,然后将结果反馈给Mathematica的 图形图。这是结果(点击放大)。
根据我要花多少时间,合理的下一步是雇用一些 降维 技术将其投影到2维并查看生成的图片(希望能提供更多信息)。但是,即使上面的简单过程也会产生一些有趣的信息:
  1. 通常,此图有助于围绕使用的某些词法转换建立直觉。
  2. 哈希标签用于执行以下操作:指示要遵循的意愿和意图进行交互。
  3. 对于英语和葡萄牙语的读者来说,星座运势显然都很重要。但是用来描述每个符号的实际单词模式却非常相似。 :)
  4. 有几种方法可以限定潜在的煽动性言论,以表明寻求真相的动机(``真实对话''集群)。
  5. 我最喜欢的#sex和#porn群集。在互联网上,它们确实是相同的。

2011年6月2日,星期四

尺寸分析和梯度下降

考虑根据损失函数$ l $,\ [
\ begin {aligned}
p ^ {(t)}&= w ^ {(t)\ top} x ^ {(t)},\\
w ^ {(t + 1)}&= w ^ {(t)}-\ eta \ frac {\ partial l} {\ partial p} x ^ {(t)}。
\ end {aligned}
\]假设最佳权重向量为$ w ^ * $。如果所有输入均按恒定量$ r $进行缩放,则最佳权重向量为$ w ^ * / r $,如果学习算法生成的输出尊重该身份,那将是很好的选择。沿着这些思路思考的更一般的方法是尺寸分析。

假设预测的维度为\ rho $,输入的维度为\ chi \。那么权重必须具有$(\ rho / \ chi)$的维数才能计算出预测方程。假设损失的大小为$ \ lambda $,这意味着学习率$ \ eta $必须具有$ \ rho ^ 2 / /(\ lambda \ chi ^ 2)$的单位才能计算出更新方程。在实践中,这意味着,如果我们有一个权重向量序列收敛到我们喜欢的某个位置,然后我们更改所有输入以使它们加倍,则学习率必须四分之一,以使生成的权重向量的整个序列将其减半,以使整个预测序列相同。

所以这些想法已经被纳入 Vowpal兔子 (实际上,这就是我意识到它们的方式:我请Vowpal团队帮助我理解我在源代码中看到的内容)。特别是在命令行上指定的$ \ eta $由$ x ^ {(t)\ top} x ^ {(t)} $规范化,对应于$ {1 / / chi ^ 2)$部分。 $ \ eta $维数;但是,这不能解决$ \ rho ^ 2 / \ lambda $部分。 (为什么这么重要?假设我创建了一个新的损失函数,该函数是另一个损失函数的输出的两倍;在命令行上指定的学习率必须降低一半。)

在处理二元模型时,我不得不弄清楚如何对此进行概括,因此我开始考虑它。通过$ x ^ {(t)\ top} x ^ {(t)} $进行归一化肯定会针对输入的全局缩放进行调整,但是如果仅对输入的一个维度进行缩放会怎样?这开始进入以下领域 预处理 这导致了自适应方法 杜契,哈桑和歌手 又名DHS(也同时在 麦克马汉和史特勒 但我将专注于DHS)。在那里他们定义了一个矩阵$ G ^ {{(t)} = \ mathrm {diag} \ left(1 / \ sqrt {\ sum_ {s \ leq t}(g ^ {{s}} _ i} ^ {2}} \ right)$,并使用更新\ [
w ^ {(t + 1)} = w ^ {(t)}-\ eta \ frac {\ partial l} {\ partial p} G ^ {(t)} x ^ {(t)},
\]其中$ g ^ {{s}} =(\ partial l / \ partial p ^ {{s}})x ^ {{s}} $。在此更新中,$ G ^ {((t)} $)的尺寸为$ \ rho /(\ lambda \ chi)\ ldots $越来越近!不幸的是,$ \ eta $仍然不是无量纲的,其尺寸为$ \ rho / \ chi $。请注意,DHS推论1中$ \ eta $的最佳选择与$ \ max_t || w_t-w ^ * || _ {\ infty} $成正比,其单位为$(\ rho / \ chi)$。换句话说,如果所有输入都加倍,我们仍然必须将以前的最佳学习率降低一半才能获得最佳行为。

这就留下了一个问题:单位$(\ rho / \ chi)$可以用来对$ \ eta $进行规范化,使得在命令行上指定的参数是无量纲的。任何与$ t $相关的变化都不在DHS的分析范围之内,但我暂时将其忽略。两件事表明了自己:$ 1 / || x ^ {(t)\ top} x ^ {(t)} || _p $和$ 1 / /(x ^ {(t)\ top} G ^ {(t)} x ^ {(t)})$。 (这些单位为$ 1 / \ chi $,但越来越近了)。它们具有不同的属性。

直观地使自适应学习算法具有优势的是,它们在频繁出现的特征上更加保守,而在很少出现的特征上更具攻击性。使用$ || x ^ {(t)\ top} x ^ {(t)} || _p $归一化,如果遇到一个例子,其中所有特征以前都被广泛地看过和训练过,那么有效学习率将很小。因此,与在训练序列中较早看到此示例相比,预测的变化将很小。相反,通过$ x ^ {(t)\ top} G ^ {(t)} x ^ {(t)} $归一化,如果遇到一个例子,其中所有特征之前都已被广泛了解和训练,则有效学习速率将被归一化以补偿,使得相对于在训练序列中较早看到该示例而言,预测的变化不会很小。另一方面,对于一个在训练序列中较晚出现的新颖特征和频繁特征混合的示例,相对于该示例是否已发生,与采用任一归一化方案的频繁特征权重相比,更新对新颖特征权重的影响将更大。在训练序列中更早。

维度$(\ rho / \ chi)$还有其他内容。自适应学习方法的关键见解之一是使用来自整个输入历史的信息,而不仅仅是当前输入,来驱动学习。到目前为止,总结所有输入的一件事是权重向量。 DHS推论1中的最佳权重向量与$ \ max_t || w_t-w ^ * || _ {\ infty} $成正比,并且由于$ w_0 = 0 $,所以这大约是$ || ww ** ||。 _ \ infty $。 $ w ^ * $的一个近似值(可能太可怕了!)是当前的权重向量。当它为零时,这在开始时是一个特别糟糕的近似值,因此我考虑将提供的$ \ eta $缩放为$ \ max(1,|| w ^ {(t)} || __infty)$,其中仅考虑$ w ^ {(t)} $的具有$ x ^ {(t)} $的相应非零值的分量。

一个实验

有很多想法,所以我决定进行一个实验。我有数以千万计的推文,这些推文根据写该推文的人的性别进行了标记。由于推文中令牌数量的变化,输入向量具有有趣的范数。我使用恒定的学习速率(vowpal按$ x ^ \ top x $进行缩放)和DHS自适应学习算法按不同的值进行缩放进行了训练。我只对数据进行了一次传递,并且报告了训练集上渐进的验证折损。我在命令行上使用$ \ eta = 1 $进行了所有这些测试。 \ [
\ begin {array} {c | c | c}
\ mbox {方法}&\ mbox {损失}&\ mbox {评论} \\ \ hline
\ mbox {最佳常数}&\ mbox {0.722}&\ mbox {在Twitter上,女人比男人多(为什么?)}} \\
\ mbox {非自适应}&\ mbox {0.651}&\ mbox {在这种情况下,大众汽车通过$ || x ^ {(t)} || ^ 2_2 $归一化}} \\
\ mbox {自适应$ || x ^ {(t)} || _1 $规范化}&\ mbox {0.588}&\\
\ mbox {自适应$ || x ^ {(t)} || _2 $规范化}&\ mbox {0.554}和\\
\ mbox {自适应$ || x ^ {(t)} || __infty $归一化}&\ mbox {0.579}&\\
\ mbox {自适应$ x ^ {(t)\ top} G ^ {(t)} x ^ {(t)} $规范化}&\ mbox {0.621}&\ mbox {比替代方案差得多! } \\
\ mbox {自适应$ \ max(1,|| w ^ {(t)} || __infty)$缩放}&\ mbox {0.579}&\\
\ mbox {自适应$ \ max(0.1,|| w ^ {(t)} || __infty)$缩放}&\ mbox {0.562}&\\
\ mbox {自适应$ \ max(1,|| w ^ {(t)} || __infty)$缩放}&\ mbox {0.579}&\ mbox {忽略$ || w ^ {{t )} || __infty $} \\
\ mbox {自适应$ \ max(0.1,|| w ^ {(t)} || __infty)$缩放}&\ mbox {0.560}&\ mbox {忽略$ || w ^ {{t )} || __infty $} \\
\ end {array}
\]
实际上,对于$ z \ leq 0.1 $尝试$ \ max(z,|| w ^ {(t)} || __infty)$会得到相同的结果。我认为这可能是由于常量特征的权重(始终以值1表示;因此,在某种意义上,它具有固定的比例,不会随输入而变化),因此我尝试在计算时不使用常量特征$ || w ^ {(t)} || _ \ infty $,但结果大致相同。

因此,该实验表明,通过自适应范数$ x ^ {(t)\ top} G ^ {(t)} x ^ {(t)} $进行归一化确实会破坏自适应策略的有效性。否则,在此实验的基础上,很难真正偏爱一种策略。话虽如此,我个人最喜欢的是$ \ max(0.1,|| w ^ {(t)} || __ \ infty)$缩放比例,因为它直接来自遗憾分析并具有正确的单位。

现在,我需要回到最初的目标:在训练一个 二元模型.