lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.6k stars 734 forks source link

Higher order derivatives by using dde.grad.jacobian and dde.grad.hessian #1580

Open svkarash opened 9 months ago

svkarash commented 9 months ago

Dear @lululxvi I have specified all essential derivatives as below. However, upon compiling model, an error occurred.

PDE is defined as

def pde2(X,U):
  x,z,t=X[:, 0:1],X[:, 1:2],X[:, 2:]
  psi,eta=U[:, 0:1], U[:, 1:]

  psi_x = dde.grad.jacobian(psi, x)
  psi_xx = dde.grad.hessian(psi, x, i=0, j=0)
  psi_xxx = dde.grad.jacobian(psi_xx, x)
  psi_xxxx = dde.grad.hessian(psi_xx, x, i=0, j=0)

  psi_z = dde.grad.jacobian(psi, z)
  psi_zz = dde.grad.hessian(psi, z, i=0, j=0)
  psi_zzz = dde.grad.jacobian(psi_zz, z)
  psi_zzzz = dde.grad.hessian(psi_zz, z, i=0, j=0)
  .....
  return ......

Then

data = dde.data.TimePDE(
    geomtime,
    pde2,
    [left_bc_psi,right_bc_psi,left_bc_eta,right_bc_eta,
     bottom_bc_psi_x, bottom_bc_psi_z,
     surface_bc1, surface_bc2, surface_bc3,
     ic_psi, ic_eta],
    num_domain=800,
    num_boundary=1000,
    num_initial=1000,
    num_test=800
)
net = dde.nn.FNN([3] + [50] * 5 + [2], "tanh", "Glorot uniform")
model.compile("adam", lr=0.001)
model.train(iterations=1000)

The Error is

AttributeError                            Traceback (most recent call last)
[<ipython-input-144-34eaae603f88>](https://localhost:8080/#) in <cell line: 2>()
      1 # Train the model
----> 2 model.compile("adam", lr=0.001)
      3 model.train(iterations=1000)
      4 # model.compile("L-BFGS-B")
      5 # losshistory, train_state = model.train(iterations=1000) #epochs

10 frames
[/usr/local/lib/python3.10/dist-packages/deepxde/utils/internal.py](https://localhost:8080/#) in wrapper(*args, **kwargs)
     20     def wrapper(*args, **kwargs):
     21         ts = timeit.default_timer()
---> 22         result = f(*args, **kwargs)
     23         te = timeit.default_timer()
     24         if config.rank == 0:

[/usr/local/lib/python3.10/dist-packages/deepxde/model.py](https://localhost:8080/#) in compile(self, optimizer, lr, loss, metrics, decay, loss_weights, external_trainable_variables)
    135 
    136         if backend_name == "tensorflow.compat.v1":
--> 137             self._compile_tensorflow_compat_v1(lr, loss_fn, decay)
    138         elif backend_name == "tensorflow":
    139             self._compile_tensorflow(lr, loss_fn, decay)

[/usr/local/lib/python3.10/dist-packages/deepxde/model.py](https://localhost:8080/#) in _compile_tensorflow_compat_v1(self, lr, loss_fn, decay)
    184             return losses
    185 
--> 186         losses_train = losses(self.data.losses_train)
    187         losses_test = losses(self.data.losses_test)
    188         total_loss = tf.math.reduce_sum(losses_train)

[/usr/local/lib/python3.10/dist-packages/deepxde/model.py](https://localhost:8080/#) in losses(losses_fn)
    170         def losses(losses_fn):
    171             # Data losses
--> 172             losses = losses_fn(
    173                 self.net.targets, self.net.outputs, loss_fn, self.net.inputs, self
    174             )

[/usr/local/lib/python3.10/dist-packages/deepxde/data/data.py](https://localhost:8080/#) in losses_train(self, targets, outputs, loss_fn, inputs, model, aux)
     11     def losses_train(self, targets, outputs, loss_fn, inputs, model, aux=None):
     12         """Return a list of losses for training dataset, i.e., constraints."""
---> 13         return self.losses(targets, outputs, loss_fn, inputs, model, aux=aux)
     14 
     15     def losses_test(self, targets, outputs, loss_fn, inputs, model, aux=None):

[/usr/local/lib/python3.10/dist-packages/deepxde/data/pde.py](https://localhost:8080/#) in losses(self, targets, outputs, loss_fn, inputs, model, aux)
    138         if self.pde is not None:
    139             if get_num_args(self.pde) == 2:
--> 140                 f = self.pde(inputs, outputs_pde)
    141             elif get_num_args(self.pde) == 3:
    142                 if self.auxiliary_var_fn is None:

[<ipython-input-128-282feef19924>](https://localhost:8080/#) in pde2(X, U)
      4 
      5   psi_x = dde.grad.jacobian(psi, x)
----> 6   psi_xx = dde.grad.hessian(psi, x)
      7   psi_xxx = dde.grad.jacobian(psi_xx, x)
      8   psi_xxxx = dde.grad.hessian(psi_xx, x)

[/usr/local/lib/python3.10/dist-packages/deepxde/gradients.py](https://localhost:8080/#) in hessian(ys, xs, component, i, j, grad_y)
    281         H[`i`][`j`].
    282     """
--> 283     return hessian._Hessians(ys, xs, component=component, i=i, j=j, grad_y=grad_y)
    284 
    285 

[/usr/local/lib/python3.10/dist-packages/deepxde/gradients.py](https://localhost:8080/#) in __call__(self, y, xs, component, i, j, grad_y)
    248             key = (id(y[0]), id(xs), component)
    249         if key not in self.Hs:
--> 250             self.Hs[key] = Hessian(y, xs, component=component, grad_y=grad_y)
    251         return self.Hs[key](i, j)
    252 

[/usr/local/lib/python3.10/dist-packages/deepxde/gradients.py](https://localhost:8080/#) in __init__(self, y, xs, component, grad_y)
    222         if grad_y is None:
    223             grad_y = jacobian(y, xs, i=component, j=None)
--> 224         self.H = Jacobian(grad_y, xs)
    225 
    226     def __call__(self, i=0, j=0):

[/usr/local/lib/python3.10/dist-packages/deepxde/gradients.py](https://localhost:8080/#) in __init__(self, ys, xs)
     20 
     21         if backend_name in ["tensorflow.compat.v1", "tensorflow", "pytorch", "paddle"]:
---> 22             self.dim_y = ys.shape[1]
     23         elif backend_name == "jax":
     24             # For backend jax, a tuple of a jax array and a callable is passed as one of

AttributeError: 'NoneType' object has no attribute 'shape'

As you can see the error refers to ----> 6 psi_xx = dde.grad.hessian(psi, x)

How can i solve the error

praksharma commented 9 months ago

Something is wrong. This is the error.

[<ipython-input-128-282feef19924>](https://localhost:8080/#) in pde2(X, U)
      4 
      5   psi_x = dde.grad.jacobian(psi, x)
----> 6   psi_xx = dde.grad.hessian(psi, x)
      7   psi_xxx = dde.grad.jacobian(psi_xx, x)
      8   psi_xxxx = dde.grad.hessian(psi_xx, x)

and this is the code you've provided.

  psi_x = dde.grad.jacobian(psi, x)
  psi_xx = dde.grad.hessian(psi, x, i=0, j=0)
  psi_xxx = dde.grad.jacobian(psi_xx, x)
  psi_xxxx = dde.grad.hessian(psi_xx, x, i=0, j=0)

Please provide the correct error. Also, it will be helpful if yiou provide the entire code.

svkarash commented 9 months ago

Dear @praksharma, you are right. However, those are the same as I read manual and command description. Assume that the command is psi_xx = dde.grad.hessian(psi, x). Since I was working on the code, some changes occurred.

praksharma commented 9 months ago

Hey I think I found the issue. According to the definition of Jacobian and Hessian here, you need to use I,j to decide which column of x and y you want to use in the differentiation.

Instead use i and j specify the column number in the tensor.

Here is a simple example,

So, don't redefine the X,U in x,z,t and psi, eta. Just use the tensor slice in the Jacobian and Hessian. For example, psi_x = dde.grad.jacobian(U,X, i=0, j=0) psi_xx = dde.grad.hessian(U, X, i=0, j=0)

Psi_xx can also be calculated as psi_xx = dde.grad.jacobian(psi_x, X, i=0, j=0) And so on.

I think slicing X and U as x,z,t and eta, psi detaches the tensor from the graph. If you don't understand this line, just ignore it.

svkarash commented 9 months ago

Hey I think I found the issue. According to the definition of Jacobian and Hessian here, you need to use I,j to decide which column of x and y you want to use in the differentiation.

Instead use i and j specify the column number in the tensor.

Here is a simple example,

  • for psi_x, i=0, j=0
  • For psi_z, i=0, j=1
  • for psi_t, i=0, j=2 Similarly,
  • for eta_x, i=1, j=0
  • For eta_z, i=1, j=1
  • for eta_t, i=1, j=2 You see i and j are column index of the X and U. This should be done for both Jacobian and Hessian.

So, don't redefine the X,U in x,z,t and psi, eta. Just use the tensor slice in the Jacobian and Hessian. For example, psi_x = dde.grad.jacobian(U,X, i=0, j=0) psi_xx = dde.grad.hessian(U, X, i=0, j=0)

Psi_xx can also be calculated as psi_xx = dde.grad.jacobian(psi_x, X, i=0, j=0) And so on.

I think slicing X and U as x,z,t and eta, psi detaches the tensor from the graph. If you don't understand this line, just ignore it.

Dear @praksharma
you are correct. I have also adjusted the definition of derivatives according to your explanation regarding X and U. It works. But this structure is not human readable :)

However, I'm encountering an issue with defining independent and dependent variables. In this scenario, all tensors have a single column, essentially forming a vector. Considering this, it seems logical that using i=0 and j=0 should be effective. This approach is also referenced in the manual for the command. I'm keen to understand why the method of defining independent and dependent variables is not proving successful.

praksharma commented 9 months ago

Any PDE is combination of derivative of output layers with respect to input layers. So if you are solving 2D fluid flow, the inputs will be x and y, whereas the output will be u,v,p. This is where you need i,j. A neural network is represented as a graph. This is how DeepXDE is designed.

If you don't like this, you can use the good old torch.autograd.grad(). You can then calculate the gradient as

psi_x_y = torch.autograd.grad(U, X, grad_outputs=torch.ones_like(U), create_graph=True)[0]
psi_x = psi_x_y[:,0]
psi_x = psi_x_y[:,1]

You see, DeepXDE abstracts all these things so you can focus on the problem rather than programming.

svkarash commented 9 months ago

Dear @praksharma , Thanks for clear explanation. Consider X=[x,z,t] and U=[psi,eta] Therefore, psi_xx = dde.grad.hessian(U, X, component=0, i=0, j=0) Now to calculate higher order derivatives such as psi_xxxx or psi_xxzz one can write

psi_xxxx = dde.grad.hessian(psi_xx, X, i=0, j=0) psi_xxzz = dde.grad.hessian(psi_xx, X, i=1, j=1)

in these cases, psi_xx is used and component=None. I believe this structure is the same as my definition in the issue. next question is how I can reshape psi_xx = dde.grad.hessian(psi, x) to be the same as psi_xx = dde.grad.hessian(U, X, component=0, i=0, j=0) Actually the error refers to the acceptable shape.

praksharma commented 9 months ago

Yes, you just need to read the first line of the docs here. It says,

H[i][j] = d^2y / dx_i dx_j

With component you can decide which y i.e, the output unit you are picking. Here are more example,

Given

psi_xx = dde.grad.hessian(U, X, component=0, i=0, j=0) 

We can write,

psi_xxxx = dde.grad.hessian(psi_xx, X, component=0, i=0, j=0)
psi_xxzz = dde.grad.hessian(psi_xx, X, component=0, i=1, j=1)
eta_xx = dde.grad.hessian(U, X, component=1, i=0, j=0)
eta_xxzt = dde.grad.hessian(eta_xx, X, component=1, i=1, j=2)

I hope this covers all the cases. Just remember in jacobian, the definition changes (here) to J[i][j] = dy_i/dx_j. So there is no component.