同时处理大数和分钟数

Question

我有一个mxn矩阵A，有非常小的或非常大的元素。例如：

A = -1e4*randn(10,20);

我想创建一个同样大小的新矩阵C，具体如下：

首先，定义一个矩阵B，其元素是A的元素的指数：

B = exp(A)

然后，定义矩阵C，使得C的每一列与B的对应列成正比，C的每一列的和等于1。换句话说，如果我们取B的一个元素除以同一列的所有元素之和，则得到C的相应元素：

C = bsxfun(@rdivide,B,sum(B));

数学上：

​

​

显然，C的所有元素都在0到1之间。但是，当使用以下代码进行计算时，所获得的矩阵包含许多NaNs：

A = -1e4*randn(10,20);

B = exp(A);

C = bsxfun(@rdivide,B,sum(B));

sum(sum(isnan(ee)))

我希望有人能提出解决这个问题的办法。先谢谢你。

更新:请注意目标是矩阵C。我定义矩阵B只是为了解释，我们可能不需要计算它。)我们其实不应该。正如@ out 89所指出的，B只包含Inf和0)。

更新2:感谢@EJG89 89链接到对数指数和技术，这可能很有用。我正在为我的问题寻找类似的分析技巧。

更新3:按照@和@EJG89 89的建议，我们可以用常数减去每一列，以便在合理的范围内获得一个新的矩阵。很明显，我们

​

对于任何常数C，我们可以选择C作为每一列的最大值：

​

(a_{max，j}是第j列的最大值)，则

​

--我觉得这个选择可以给出一个很好的近似，但我不知道它有多好？

新的守则是：

A = bsxfun(@minus,A,max(A));
B = exp(A);
C = bsxfun(@rdivide,B,sum(B));

Answer 1

直接的方法是不可能的。

A = -1e4*randn(10,20);
B = exp(A);

问题是B已经产生Inf或0。

屈服除法不可能还南的。

ndoublemax = realmax
nsinglemax = realmax('single')

ndoublemin = realmin
nsinglemin = realmin('single')

屈服：

ndoublemax =

1.7977e+308

nsinglemax =

3.4028e+038

ndoublemin =

2.2251e-308

nsinglemin =

1.1755e-038

所以你的B值明显超过了那些值。

一个近似可以通过用：http://lingpipe-blog.com/2009/06/25/log-sum-of-exponentials/重写指数函数来实现。

所以首先重写: log(Cij) = log (exp(Aij)/sumk(exp(A(Kj)

log(Cij) = Aij - log( sumk(exp(A(Kj)

其中最后一项可以用:m= maxk(A(kj))来逼近。

log(sumk(exp(A(Kj)=m+ log(exp(A(kj)-m))

这使得MatLab能够处理较小的值，也就是说，规范化。但是，仍然没有确切的答案，因为一些值变为零，消除了它们有时很重要的贡献。

致以亲切的问候，

恩斯特·扬

Answer 2

在讨论@ all 89之后，您可以用最大可表示数替换无限值，然后将eps添加到分母中的所有数字中。这样，就可以避免Inf。

A = -1e4*randn(10,20);

B = exp(A);
B(isinf(B)) = realmax;

C = bsxfun(@rdivide,B,sum(B)+eps);

sum(isnan(C(:)))

Answer 3

一种方法可以是诉诸符号计算：

A = -1e4*randn(10,20);
S = sym(A,'d'); %// convert to symbolic
C = exp(S)./repmat(sum(exp(S)),size(S,1),1);

以符号形式给出了结果C。例如，第一列将类似于

>> C(:,1)
ans =
   2.9347694526167482187885232843876*10^(-790)
  1.548257079197982942994953632259*10^(-10497)
  8.4257350773369612716806306890481*10^(-6390)
  3.1937456016712092927990884814437*10^(-7937)
  4.0377045095817442881784227420794*10^(-4076)
                                           1.0
  3.3704875581881026436375125643672*10^(-6482)
 8.9221015470879963245101895439354*10^(-12246)
 9.4968643813486260650483647051531*10^(-11252)
 1.8777491443104052522832960625121*10^(-11084)

但是请注意，即使用符号计算，最大的项也是1，而实际上它(稍微)更小。当然，如果然后转换为双倍，则只能得到0或1的值：

>> Cd = eval(C);
>> Cd(:,1)
ans =
     0
     0
     0
     0
     0
     1
     0
     0
     0
     0