标签: hessian-matrix

我想查找给定图像的山脊。（边缘没有边缘！）示例如下图所示

我认为Hessian矩阵会直观地工作。因此，我从2D-高斯方程开始对Hessian矩阵内核进行了硬编码，如下所述。 如何构建2D粗麻布矩阵内核

我使用surf可视化的方法创建了3个二阶导数内核（D_xx，D_yy和D_xy），它们看上去都是正确的。

然后，我应用了这些内核并对图像进行了2D卷积。

我不确定下一步该怎么做，是否需要使用D_xx，D_yy和D_xy来表示特征值和向量？我们如何通过对每个像素进行2×2矩阵的特征分析来从图像中提取出棱线？任何想法，公式甚至代码都将大有帮助。

附带的代码生成二维Hessian矩阵

[x y]=meshgrid(round(-N/2):round(N/2), round(-N/2):round(N/2));
common = x.^2+y.^2;
Lxx = ((-1+x.^2/(sigma^2)).*exp(-common/(2*sigma^2))) / (2*pi*sigma^4);
Lxx = Lxx./ sum(Lxx(:));

Lyy = ((-1+y.^2/(sigma^2)).*exp(-common/(2*sigma^2))) / (2*pi*sigma^4);
Lyy = Lyy./ sum(Lyy(:));

Lxy = ((x.*y)/(2*pi*sigma^6)).*exp(-common/(2*sigma^2));
Lxy = Lxy./ sum(Lxy(:));

我正在尝试使用Theano来计算关于向量以及几个标量的函数的粗糙度(编辑:也就是说,我基本上希望附加到我正在计算粗麻布的向量的标量) .这是一个最小的例子:

import theano
import theano.tensor as T
A = T.vector('A')
b,c = T.scalars('b','c')
y = T.sum(A)*b*c

我的第一次尝试是:

hy = T.hessian(y,[A,b,c])

哪个失败了 AssertionError: tensor.hessian expects a (list of) 1 dimensional variable as 'wrt'

我的第二次尝试是将A,b和c与:

wrt = T.concatenate([A,T.stack(b,c)])
hy = T.hessian(y,[wrt])

哪个失败了 DisconnectedInputError: grad method was asked to compute the gradient with respect to a variable that is not part of the computational graph of the cost, or is used only by a non-differentiable operator: Join.0

在这种情况下计算粗麻线的正确方法是什么？

更新:为了澄清我在寻找什么,假设A是2元素向量.然后黑森州将是:

[[d2y/d2A1, d2y/dA1dA2, d2y/dA1dB, …

我正在尝试获得 theano 梯度和 hessian 的矢量化版本，即我想在几个点计算梯度和 hessian，在矩阵中给出，如下所示：

我有一个功能：

f(x_1,x_2,..,x_n)=exp(x_1^2+x_2^2+...+x_n^2)

我想用一个命令在多个点计算它的梯度。我可以这样做：

x = T.matrix('x')
y = T.diag(T.exp(T.dot(x,x.T)))
J = theano.grad(cost = y.sum(), wrt = x)
f = theano.function(inputs = [x], outputs = J)
f([[1,2],[3,4]])

它返回一个矩阵，其中的行是在点 (1,2) 和 (3,4) 处计算的梯度。我想为 hessian 获得相同的结果（在这种情况下，它将是一个与矩阵相反的 3 维张量，但想法相同）。以下代码：

H = theano.gradient.hessian(cost = y.sum(), wrt = x)

返回错误：

AssertionError: tensor.hessian expects a (list of) 1 dimensional variable as `wrt`

我能够使用以下代码获得适当的结果

J = theano.grad(cost = y.sum(), wrt = x)
H = theano.gradient.jacobian(expression = J.flatten(), wrt = x)
g = …

我之前曾问过这个问题,并希望继续跟进,因为我尝试了其他一些事情并且他们没有完全解决问题.

我本质上是在尝试优化R中的NLP类型问题,它具有二进制和整数约束.相同的代码如下:

# Input Data
DTM <- sample(1:30,10,replace=T)
DIM <- rep(30,10)
Price <- 100 - seq(0.4,1,length.out=10)

# Variables that shall be changed to find optimal solution
Hike <- c(1,0,0,1,0,0,0,0,0,1)
Position <- c(0,1,-2,1,0,0,0,0,0,0)

# Bounds for Hikes/Positions
HikeLB <- rep(0,10)
HikeUB <- rep(1,10)
PositionLB <- rep(-2,10)
PositionUB <- rep(2,10)

library(Rsolnp)

# x <- c(Hike, Position)
# Combining two arrays into one since I want 
# to optimize using both these variables

opt_func <- function(x) {

  Hike <- head(x,length(x)/2)
  Position …

我在Keras和TensorFlow一起工作.我有一个深度神经模型(预测自动编码器).我正在做一些与此类似的事情:https://arxiv.org/abs/1612.00796 - 我试图了解输出中给定图层中变量的影响.

为此,我需要找到关于特定层输出的损失(L)的二阶导数(Hessian):

对角线条目就足够了.L是标量,s是1乘n.

我先尝试了什么:

dLds = tf.gradients(L, s)  # works fine to get first order derivatives
d2Lds2 = tf.gradients(dLds, s)  # throws an error
TypeError: Second-order gradient for while loops not supported.

我也尝试过:

d2Lds2 = tf.hessians(L, s)
ValueError: Computing hessians is currently only supported for one-dimensional tensors. Element number 0 of `xs` has 2 dimensions.

我无法改变s的形状,因为它是神经网络的一部分(LSTM的状态).第一个维度(batch_size)已经设置为1,我认为我不能摆脱它.

我无法重塑s因为它打破了渐变的流动,例如:

tf.gradients(L, tf.reduce_sum(s, axis=0))

得到:

[None]

关于在这种情况下我该怎么做的任何想法？

library(nlme)
mydat <- 
  structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L, 
                           92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 
                           92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 
                           94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 
                           95L, 95L, 95L), days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9, 
                                                    2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9, 
                                                    17.3, 17.7, 19.4, 2.3, …

我正在努力解决以下问题。简而言之：两个不同的软件包（Aptech 的 Gauss 和 R）在最大似然法中产生完全不同的 Hessian 矩阵。我使用相同的程序（BFGS），完全相同的数据，相同的最大似然公式（这是一个非常简单的 logit 模型），具有完全相同的起始值，令人困惑的是，我得到了相同的参数和 log-结果可能性。两个程序中只有 Hessian 矩阵不同，因此标准误差的估计和统计推断不同。

在这个特定的例子中，它并没有出现太大的偏差，但是模型的每一个增加的复杂性都会增加差异，所以如果我尝试估计我的最终模型，两个程序都会产生完全错误的结果。

有谁知道，这两个程序在计算 Hessian 矩阵的方式上有何不同，以及获得相同结果的正确方法可能是什么？

编辑：在 R （高斯）代码中，向量X ( alt ) 是自变量，由两列向量组成，第一列完全是 1，第二列是受试者的响应。向量y ( itn ) 是因变量，由一列和受试者的响应组成。该示例（R 代码和数据集）取自http://www.polsci.ucsb.edu/faculty/glasgow/ps206/ps206.html，只是作为重现和隔离问题的示例。

我已附上代码（高斯和 R 语法）和输出。

任何帮助将不胜感激。谢谢 ：）

高斯：

start={ 0.95568840 , -0.20459156 };

library maxlik,pgraph;
maxset;
_max_Algorithm = 2;
_max_Diagnostic = 1;
{betaa,f,g,cov,ret} = maxlik(XMAT,0,&ll,start);
call maxprt(betaa,f,g,cov,ret);
print _max_FinalHess;

proc ll(b,XMAT);
local exb, probo, logexb, yn, logexbn, yt, ynt, logl;

    exb = EXP(alt*b);
    //print exb;
    probo = exb./(1+exb); …

因此，我对 Python 中的 ML/AI 游戏相对较新，目前正在研究围绕 XGBoost 自定义目标函数实现的问题。

我的微分方程知识相当生疏，所以我创建了一个带有梯度和 hessian 的自定义 obj 函数，该函数对均方误差函数进行建模，该函数作为 XGBRegressor 中的默认目标函数运行，以确保我正确执行所有这些操作。问题是，模型的结果（错误输出很接近，但在大多数情况下并不相同（并且在某些点上相差很大）。我不知道我做错了什么，也不知道如果我做错了什么，这怎么可能我的计算是正确的。如果你们都可以看看这个，也许可以深入了解我错在哪里，那就太棒了！

没有自定义函数的原始代码是：

    import xgboost as xgb

    reg = xgb.XGBRegressor(n_estimators=150, 
                   max_depth=2,
                   objective ="reg:squarederror", 
                   n_jobs=-1)

    reg.fit(X_train, y_train)

    y_pred_test = reg.predict(X_test)

我的 MSE 自定义目标函数如下：

    def gradient_se(y_true, y_pred):
        #Compute the gradient squared error.
        return (-2 * y_true) + (2 * y_pred)

    def hessian_se(y_true, y_pred):
        #Compute the hessian for squared error
        return 0*(y_true + y_pred) + 2

   def custom_se(y_true, y_pred):
        #squared error objective. A simplified version of MSE used as
        #objective function. …

我期待看到如何编译计算对数似然的Hessian的函数,以便它可以有效地与不同的参数集一起使用.

这是一个例子.

假设我们有一个函数来计算logit模型的对数似然,其中y是向量,x是矩阵.beta是参数的向量.

 pLike[y_, x_, beta_] :=
 Module[
  {xbeta, logDen},
  xbeta = x.beta;
  logDen = Log[1.0 + Exp[xbeta]];
  Total[y*xbeta - logDen]
  ]

鉴于以下数据,我们可以如下使用它

In[1]:= beta = {0.5, -1.0, 1.0};

In[2]:= xmat = 
  Table[Flatten[{1, 
     RandomVariate[NormalDistribution[0.0, 1.0], {2}]}], {500}];

In[3]:= xbeta = xmat.beta;

In[4]:= prob = Exp[xbeta]/(1.0 + Exp[xbeta]);

In[5]:= y = Map[RandomVariate[BernoulliDistribution[#]] &, prob] ;

In[6]:= Tally[y]

Out[6]= {{1, 313}, {0, 187}}

In[9]:= pLike[y, xmat, beta]

Out[9]= -272.721

我们可以写下它的粗麻布如下

 hessian[y_, x_, z_] :=
  Module[{},
   D[pLike[y, x, z], {z, 2}]
  ]


In[10]:= …

我在使用R中的函数arma {tseries}和arima {stats}拟合ARMA模型时发现了一些奇怪的东西.

两个函数采用的估计程序存在根本区别,即在arma {stats}中的卡尔曼滤波器,而不是arma {tseries}中的ML估计.

鉴于两个函数之间的估计过程存在差异,如果我们使用相同的时间序列,则不会期望两个函数的结果完全不同.

好吧他们可以!

生成以下时间序列并添加2个异常值.

set.seed(1010)
ts.sim <- abs(arima.sim(list(order = c(1,0,0), ar = 0.7), n = 50))
ts.sim[8] <- ts.sim[12]*8
ts.sim[35] <- ts.sim[32]*8

使用这两个函数拟合ARMA模型.

# Works perfectly fine
arima(ts.sim, order = c(1,0,0))
# Works perfectly fine
arma(ts.sim, order = c(1,0))

将时间序列的级别更改为10亿

# Introduce a multiplicative shift
ts.sim.1 <- ts.sim*1000000000
options(scipen = 999)
summary(ts.sim.1)

使用2个函数拟合ARMA模型:

# Works perfectly fine
arma(ts.sim.1, order = c(1,0))

# Does not work
arima(ts.sim.1, order = c(1,0,0))

## Error …