My answer for my question: yes, it can be shown that gradient for logistic loss is equal to difference between true values and predicted probabilities. Brief explanation was found here.

First, logistic loss is just negative log-likelihood, so we can start with expression for log-likelihood (p. 74 - this expression is log-likelihood itself, not negative log-likelihood):

is logistic function: , where is predicted values before logistic transformation (i.e., log-odds):

First derivative obtained using Wolfram Alpha:

After multiplying by :

After changing sign we have expression for gradient of logistic loss function:

The log loss function is given as:

where

Taking the partial derivative we get the gradient as

Thus we get the negative of gradient as p-y.

Similar calculations can be done to obtain the hessian.

Note that this algorithm is memory bandwidth limited. If you optimize this in a larger context (real application) there is very likely a higher speedup possible.

Example

import numpy as np

#https://stackoverflow.com/a/29863846/4045774
def sigmoid(x):  
    return np.exp(-np.logaddexp(0, -x))

def calc_loss_grad_1(weights, X_batch, y_batch):
    n_samples, n_features = X_batch.shape
    loss_grad = np.zeros((n_features,))
    for i in range(n_samples):
        sigm = sigmoid(X_batch[i,:] @ weights)
        loss_grad += - (y_batch[i] - sigm) * X_batch[i]            
    return loss_grad

def calc_loss_grad_2(weights, X_batch, y_batch):
    sigm =-y_batch+sigmoid(X_batch@weights)          
    return sigm@X_batch

weights=np.random.rand(n_features)
X_batch=np.random.rand(n_samples,n_features)
y_batch=np.random.rand(n_samples)

#n_samples=int(1e5)
#n_features=int(1e4)
%timeit res=calc_loss_grad_1(weights, X_batch, y_batch)
#1.79 s ± 35.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res=calc_loss_grad_2(weights, X_batch, y_batch)
#539 ms ± 21.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

#n_samples=int(1e3)
#n_features=int(1e2)
%timeit res=calc_loss_grad_1(weights, X_batch, y_batch)
#3.68 ms ± 44.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit res=calc_loss_grad_2(weights, X_batch, y_batch)
#49.1 µs ± 503 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

There are some small mistakes like you should use np.sum(Y*np.log(A) + (1-Y)*np.log(1-A)) / m in place of using .mean() and the next mistake that I think is replace np.subtract(A-Y) with simple A-Y bcz. there is no need for numpy in this. It's working for me.

def propagate(w, b, X, Y):
"""
Implement the cost function and its gradient for the propagation explained above

Arguments:
w -- weights, a numpy array of size (num_px * num_px * 3, 1)
b -- bias, a scalar
X -- data of size (num_px * num_px * 3, number of examples)
Y -- true "label" vector (containing 0 if non-cat, 1 if cat) of size (1, number of examples)

Return:
cost -- negative log-likelihood cost for logistic regression
dw -- gradient of the loss with respect to w, thus same shape as w
db -- gradient of the loss with respect to b, thus same shape as b

Tips:
- Write your code step by step for the propagation. np.log(), np.dot()
"""

m = X.shape[1]

# FORWARD PROPAGATION (FROM X TO COST)
### START CODE HERE ### (≈ 2 lines of code)
A = sigmoid(np.dot(w.T,X)+b)                                   # compute activation
cost = -np.sum(Y*np.log(A) + (1-Y)*np.log(1-A)) / m                              # compute cost
### END CODE HERE ###  

# BACKWARD PROPAGATION (TO FIND GRAD)
### START CODE HERE ### (≈ 2 lines of code)
dw = np.dot(X,(A-Y).T)/m
db = np.sum(A-Y,axis=1)/m
### END CODE HERE ###

assert(dw.shape == w.shape)
assert(db.dtype == float)
cost = np.squeeze(cost)
assert(cost.shape == ())

grads = {"dw": dw,
         "db": db}

return grads, cost

It's not advisable to remove the log term simply because it's monotone since intuitively as a score function the log term in policy gradient theorem (PGT) is not arbitrary but rather a necessary step to ensure proper scaling and direction of the gradients especially for actions with small probabilities. We can have a simple 1-d parametere space example to show this.

Let's say we have two possible actions $a_1,a_2$ of a simple bandit problem with the parameterized policy for both actions defined as $\pi_{\theta}(a_1)=\sigma(\theta), \pi_{\theta}(a_2)=1-\sigma(\theta)$, where $\sigma$ is the standard sigmoid function. Also suppose the rewards for the actions are $r(a_1)=1, r(a_2)=0$, since it's a bandit without state the return $G_t$ is simply $r(a)$ for each action $a$. Then it's straightforward to compute gradient as follows: $$\nabla_{\theta} J(\theta)=\mathbb{E}_{\pi_{\theta}}[\nabla_{\theta}\log\pi_{\theta}(a)G_t]=\pi_{\theta}(a_1)(1-\sigma(\theta))r(a_1)+\pi_{\theta}(a_2)(-\sigma(\theta))r(a_2)=\sigma(\theta)(1-\sigma(\theta))$$ Similarly you can get the gradient for the case of removing the log as $$\nabla_{\theta} J'(\theta)=\sigma(\theta)^2(1-\sigma(\theta))$$

Therefore clearly even in 1-d bandit case the gradient without log will become extremely small if $\sigma(\theta)$, i.e., $\pi_{\theta}(a_1)$ is very small, which means the initially very unlikely actions in the to-be-optimized policy cannot get updated efficiently consistent with above intuition due to losing the proper scaling effect of the log. And in the usual multi-dimensional policy for MDP cases both the magnitude and direction of the gradient will be similarly negatively impacted.

Finally with the common log gradient trick the existing form in PGT can be simply transformed to a quotient where the numerator is just your gradient of the policy without log and the denominator is the same policy. From the balance need of exploitation-exploration, this denominator is required for exploration since if the probability of taking certain action in a state is small, then the algorithm would update the parameters so that the probability of taking that action can increase. The other term $G_t \approx Q_t(s_t,a_t)$ reflects that if an action value is great then the algorithm intends to update the parameters so that the probability of taking that action is enhanced which is exploitation. Therefore if you remove the denominator policy then equivalently you're not exploring enough.