We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
F.mean_absolute_error numerically unstable with float16 arrays
In #5053, float16 support has been enabled for [F.mean_absolute_error](https://github.com/chainer/chainer/blob/master/chainer/functions/loss/mean_absolute_error.py), but it seems to produce NaN values quite easily. Usually this happens when using big batch sizes and/or if the absolute error difference is large.

The calculation is done by summing over all the absolute differences, and then dividing by the number of elements in the array. However, it appears that the summing can produce large numbers outside the possible range for `float16`. The current implementation can be broken down as follows:

```python
def forward_cpu(self, inputs):
    x0, x1 = inputs
    self.diff = x0 - x1
    diff = self.diff.ravel()
    abs_diff = abs(diff)
    summed_abs_diff = abs_diff.sum()  # numerically unstable, can result in inf
    mean_abs_error = np.array(summed_abs_diff / diff.size, dtype=diff.dtype)
    return mean_abs_error
```

Code to reproduce error:

```python
import chainer.functions as F
import numpy as np

a = np.full(shape=(64,1,16,16), fill_value=2, dtype=np.float16)
b = np.full(shape=(64,1,16,16), fill_value=-2, dtype=np.float16)

loss = F.mean_absolute_error(a,b)
# /home/user/.local/share/virtualenvs/.../lib/python3.6/site-packages/numpy/core/_methods.py:36: RuntimeWarning: overflow encountered in reduce
#  return umr_sum(a, axis, dtype, out, keepdims, initial)
# variable(inf)

loss = F.mean_absolute_error(a.astype("float32"), b.astype("float32"))
# variable(4.)
```

Note that the actual loss (4) would still be valid in the float16 range, it is just that summing over many values results in an `inf`, which cannot then be divided to get a proper number.

Workaround ideas:

I've noticed the new `mixed16` mode that was implemented in #6456, and was wondering if there might be a similar way to do the intermediate calculations in `float32`, and cast the result back into `float16`? Thoughts?

System info:
```
Platform: Linux-4.15.0-46-generic-x86_64-with-debian-buster-sid
Chainer: 6.0.0b3
NumPy: 1.16.2
CuPy:
  CuPy Version          : 6.0.0b3
  CUDA Root             : /usr/local/cuda
  CUDA Build Version    : 10000
  CUDA Driver Version   : 10010
  CUDA Runtime Version  : 10000
  cuDNN Build Version   : 7402
  cuDNN Version         : 7402
  NCCL Build Version    : 2307
  NCCL Runtime Version  : 2307
iDeep: Not Available
```

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `chainer/functions/loss/mean_absolute_error.py`
Content:
```
1 import numpy
2 
3 import chainer
4 from chainer import backend
5 from chainer import function_node
6 from chainer.utils import type_check
7 
8 
9 class MeanAbsoluteError(function_node.FunctionNode):
10 
11     """Mean absolute error function."""
12 
13     def check_type_forward(self, in_types):
14         type_check._argname(in_types, ('x0', 'x1'))
15         type_check.expect(
16             in_types[0].dtype.kind == 'f',
17             in_types[0].dtype == in_types[1].dtype,
18             in_types[0].shape == in_types[1].shape
19         )
20 
21     def forward_cpu(self, inputs):
22         x0, x1 = inputs
23         self.diff = x0 - x1
24         diff = self.diff.ravel()
25         return numpy.array(abs(diff).sum() / diff.size, dtype=diff.dtype),
26 
27     def forward_gpu(self, inputs):
28         x0, x1 = inputs
29         self.diff = x0 - x1
30         diff = self.diff.ravel()
31         return abs(diff).sum() / diff.dtype.type(diff.size),
32 
33     def backward(self, indexes, grad_outputs):
34         gy, = grad_outputs
35         coeff = gy * gy.data.dtype.type(1. / self.diff.size)
36         coeff = chainer.functions.broadcast_to(coeff, self.diff.shape)
37         gx0 = coeff * backend.get_array_module(gy.data).sign(self.diff)
38         return gx0, -gx0
39 
40 
41 def mean_absolute_error(x0, x1):
42     """Mean absolute error function.
43 
44     The function computes the mean absolute error between two variables. The
45     mean is taken over the minibatch. Args ``x0`` and ``x1`` must have the
46     same dimensions. This function first calculates the absolute value
47     differences between the corresponding elements in x0 and x1, and then
48     returns the mean of those differences.
49 
50     Args:
51         x0 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.
52         x1 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.
53 
54     Returns:
55         ~chainer.Variable:
56             A variable holding an array representing the mean absolute
57             error of two inputs.
58 
59     .. admonition:: Example
60 
61         1D array examples:
62 
63         >>> x = np.array([1, 2, 3]).astype(np.float32)
64         >>> y = np.array([0, 0, 0]).astype(np.float32)
65         >>> F.mean_absolute_error(x, y)
66         variable(2.)
67         >>> x = np.array([1, 2, 3, 4, 5, 6]).astype(np.float32)
68         >>> y = np.array([7, 8, 9, 10, 11, 12]).astype(np.float32)
69         >>> F.mean_absolute_error(x, y)
70         variable(6.)
71 
72         2D array example:
73 
74         In this example, there are 4 elements, and thus 4 errors
75         >>> x = np.array([[1, 2], [3, 4]]).astype(np.float32)
76         >>> y = np.array([[8, 8], [8, 8]]).astype(np.float32)
77         >>> F.mean_absolute_error(x, y)
78         variable(5.5)
79 
80         3D array example:
81 
82         In this example, there are 8 elements, and thus 8 errors
83         >>> x = np.reshape(np.array([1, 2, 3, 4, 5, 6, 7, 8]), (2, 2, 2))
84         >>> y = np.reshape(np.array([8, 8, 8, 8, 8, 8, 8, 8]), (2, 2, 2))
85         >>> x = x.astype(np.float32)
86         >>> y = y.astype(np.float32)
87         >>> F.mean_absolute_error(x, y)
88         variable(3.5)
89 
90     """
91     return MeanAbsoluteError().apply((x0, x1))[0]
92
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/chainer/functions/loss/mean_absolute_error.py b/chainer/functions/loss/mean_absolute_error.py
--- a/chainer/functions/loss/mean_absolute_error.py
+++ b/chainer/functions/loss/mean_absolute_error.py
@@ -6,6 +6,15 @@
 from chainer.utils import type_check
 
 
+def _get_intermediate_dtype(dtype):
+    # Returns the dtype for intermediate calculation.
+    # For float16 input, float32 is used.
+    # Otherwise the same dtype as the parameter is used.
+    if dtype == numpy.float16:
+        return numpy.float32
+    return dtype
+
+
 class MeanAbsoluteError(function_node.FunctionNode):
 
     """Mean absolute error function."""
@@ -21,14 +30,19 @@
     def forward_cpu(self, inputs):
         x0, x1 = inputs
         self.diff = x0 - x1
-        diff = self.diff.ravel()
-        return numpy.array(abs(diff).sum() / diff.size, dtype=diff.dtype),
+        orig_dtype = self.diff.dtype
+        dtype = _get_intermediate_dtype(orig_dtype)
+        diff = self.diff.ravel().astype(dtype, copy=False)
+        return numpy.array(abs(diff).sum() / diff.size, dtype=orig_dtype),
 
     def forward_gpu(self, inputs):
         x0, x1 = inputs
         self.diff = x0 - x1
-        diff = self.diff.ravel()
-        return abs(diff).sum() / diff.dtype.type(diff.size),
+        orig_dtype = self.diff.dtype
+        dtype = _get_intermediate_dtype(orig_dtype)
+        diff = self.diff.ravel().astype(dtype, copy=False)
+        return (abs(diff).sum() / diff.dtype.type(diff.size)).astype(
+            orig_dtype, copy=False),
 
     def backward(self, indexes, grad_outputs):
         gy, = grad_outputs

{"golden_diff": "diff --git a/chainer/functions/loss/mean_absolute_error.py b/chainer/functions/loss/mean_absolute_error.py\n--- a/chainer/functions/loss/mean_absolute_error.py\n+++ b/chainer/functions/loss/mean_absolute_error.py\n@@ -6,6 +6,15 @@\n from chainer.utils import type_check\n \n \n+def _get_intermediate_dtype(dtype):\n+    # Returns the dtype for intermediate calculation.\n+    # For float16 input, float32 is used.\n+    # Otherwise the same dtype as the parameter is used.\n+    if dtype == numpy.float16:\n+        return numpy.float32\n+    return dtype\n+\n+\n class MeanAbsoluteError(function_node.FunctionNode):\n \n     \"\"\"Mean absolute error function.\"\"\"\n@@ -21,14 +30,19 @@\n     def forward_cpu(self, inputs):\n         x0, x1 = inputs\n         self.diff = x0 - x1\n-        diff = self.diff.ravel()\n-        return numpy.array(abs(diff).sum() / diff.size, dtype=diff.dtype),\n+        orig_dtype = self.diff.dtype\n+        dtype = _get_intermediate_dtype(orig_dtype)\n+        diff = self.diff.ravel().astype(dtype, copy=False)\n+        return numpy.array(abs(diff).sum() / diff.size, dtype=orig_dtype),\n \n     def forward_gpu(self, inputs):\n         x0, x1 = inputs\n         self.diff = x0 - x1\n-        diff = self.diff.ravel()\n-        return abs(diff).sum() / diff.dtype.type(diff.size),\n+        orig_dtype = self.diff.dtype\n+        dtype = _get_intermediate_dtype(orig_dtype)\n+        diff = self.diff.ravel().astype(dtype, copy=False)\n+        return (abs(diff).sum() / diff.dtype.type(diff.size)).astype(\n+            orig_dtype, copy=False),\n \n     def backward(self, indexes, grad_outputs):\n         gy, = grad_outputs\n", "issue": "F.mean_absolute_error numerically unstable with float16 arrays\nIn #5053, float16 support has been enabled for [F.mean_absolute_error](https://github.com/chainer/chainer/blob/master/chainer/functions/loss/mean_absolute_error.py), but it seems to produce NaN values quite easily. Usually this happens when using big batch sizes and/or if the absolute error difference is large.\r\n\r\nThe calculation is done by summing over all the absolute differences, and then dividing by the number of elements in the array. However, it appears that the summing can produce large numbers outside the possible range for `float16`. The current implementation can be broken down as follows:\r\n\r\n```python\r\ndef forward_cpu(self, inputs):\r\n    x0, x1 = inputs\r\n    self.diff = x0 - x1\r\n    diff = self.diff.ravel()\r\n    abs_diff = abs(diff)\r\n    summed_abs_diff = abs_diff.sum()  # numerically unstable, can result in inf\r\n    mean_abs_error = np.array(summed_abs_diff / diff.size, dtype=diff.dtype)\r\n    return mean_abs_error\r\n```\r\n\r\nCode to reproduce error:\r\n\r\n```python\r\nimport chainer.functions as F\r\nimport numpy as np\r\n\r\na = np.full(shape=(64,1,16,16), fill_value=2, dtype=np.float16)\r\nb = np.full(shape=(64,1,16,16), fill_value=-2, dtype=np.float16)\r\n\r\nloss = F.mean_absolute_error(a,b)\r\n# /home/user/.local/share/virtualenvs/.../lib/python3.6/site-packages/numpy/core/_methods.py:36: RuntimeWarning: overflow encountered in reduce\r\n#  return umr_sum(a, axis, dtype, out, keepdims, initial)\r\n# variable(inf)\r\n\r\nloss = F.mean_absolute_error(a.astype(\"float32\"), b.astype(\"float32\"))\r\n# variable(4.)\r\n```\r\n\r\nNote that the actual loss (4) would still be valid in the float16 range, it is just that summing over many values results in an `inf`, which cannot then be divided to get a proper number.\r\n\r\nWorkaround ideas:\r\n\r\nI've noticed the new `mixed16` mode that was implemented in #6456, and was wondering if there might be a similar way to do the intermediate calculations in `float32`, and cast the result back into `float16`? Thoughts?\r\n\r\nSystem info:\r\n```\r\nPlatform: Linux-4.15.0-46-generic-x86_64-with-debian-buster-sid\r\nChainer: 6.0.0b3\r\nNumPy: 1.16.2\r\nCuPy:\r\n  CuPy Version          : 6.0.0b3\r\n  CUDA Root             : /usr/local/cuda\r\n  CUDA Build Version    : 10000\r\n  CUDA Driver Version   : 10010\r\n  CUDA Runtime Version  : 10000\r\n  cuDNN Build Version   : 7402\r\n  cuDNN Version         : 7402\r\n  NCCL Build Version    : 2307\r\n  NCCL Runtime Version  : 2307\r\niDeep: Not Available\r\n```\n", "before_files": [{"content": "import numpy\n\nimport chainer\nfrom chainer import backend\nfrom chainer import function_node\nfrom chainer.utils import type_check\n\n\nclass MeanAbsoluteError(function_node.FunctionNode):\n\n    \"\"\"Mean absolute error function.\"\"\"\n\n    def check_type_forward(self, in_types):\n        type_check._argname(in_types, ('x0', 'x1'))\n        type_check.expect(\n            in_types[0].dtype.kind == 'f',\n            in_types[0].dtype == in_types[1].dtype,\n            in_types[0].shape == in_types[1].shape\n        )\n\n    def forward_cpu(self, inputs):\n        x0, x1 = inputs\n        self.diff = x0 - x1\n        diff = self.diff.ravel()\n        return numpy.array(abs(diff).sum() / diff.size, dtype=diff.dtype),\n\n    def forward_gpu(self, inputs):\n        x0, x1 = inputs\n        self.diff = x0 - x1\n        diff = self.diff.ravel()\n        return abs(diff).sum() / diff.dtype.type(diff.size),\n\n    def backward(self, indexes, grad_outputs):\n        gy, = grad_outputs\n        coeff = gy * gy.data.dtype.type(1. / self.diff.size)\n        coeff = chainer.functions.broadcast_to(coeff, self.diff.shape)\n        gx0 = coeff * backend.get_array_module(gy.data).sign(self.diff)\n        return gx0, -gx0\n\n\ndef mean_absolute_error(x0, x1):\n    \"\"\"Mean absolute error function.\n\n    The function computes the mean absolute error between two variables. The\n    mean is taken over the minibatch. Args ``x0`` and ``x1`` must have the\n    same dimensions. This function first calculates the absolute value\n    differences between the corresponding elements in x0 and x1, and then\n    returns the mean of those differences.\n\n    Args:\n        x0 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.\n        x1 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.\n\n    Returns:\n        ~chainer.Variable:\n            A variable holding an array representing the mean absolute\n            error of two inputs.\n\n    .. admonition:: Example\n\n        1D array examples:\n\n        >>> x = np.array([1, 2, 3]).astype(np.float32)\n        >>> y = np.array([0, 0, 0]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(2.)\n        >>> x = np.array([1, 2, 3, 4, 5, 6]).astype(np.float32)\n        >>> y = np.array([7, 8, 9, 10, 11, 12]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(6.)\n\n        2D array example:\n\n        In this example, there are 4 elements, and thus 4 errors\n        >>> x = np.array([[1, 2], [3, 4]]).astype(np.float32)\n        >>> y = np.array([[8, 8], [8, 8]]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(5.5)\n\n        3D array example:\n\n        In this example, there are 8 elements, and thus 8 errors\n        >>> x = np.reshape(np.array([1, 2, 3, 4, 5, 6, 7, 8]), (2, 2, 2))\n        >>> y = np.reshape(np.array([8, 8, 8, 8, 8, 8, 8, 8]), (2, 2, 2))\n        >>> x = x.astype(np.float32)\n        >>> y = y.astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(3.5)\n\n    \"\"\"\n    return MeanAbsoluteError().apply((x0, x1))[0]\n", "path": "chainer/functions/loss/mean_absolute_error.py"}], "after_files": [{"content": "import numpy\n\nimport chainer\nfrom chainer import backend\nfrom chainer import function_node\nfrom chainer.utils import type_check\n\n\ndef _get_intermediate_dtype(dtype):\n    # Returns the dtype for intermediate calculation.\n    # For float16 input, float32 is used.\n    # Otherwise the same dtype as the parameter is used.\n    if dtype == numpy.float16:\n        return numpy.float32\n    return dtype\n\n\nclass MeanAbsoluteError(function_node.FunctionNode):\n\n    \"\"\"Mean absolute error function.\"\"\"\n\n    def check_type_forward(self, in_types):\n        type_check._argname(in_types, ('x0', 'x1'))\n        type_check.expect(\n            in_types[0].dtype.kind == 'f',\n            in_types[0].dtype == in_types[1].dtype,\n            in_types[0].shape == in_types[1].shape\n        )\n\n    def forward_cpu(self, inputs):\n        x0, x1 = inputs\n        self.diff = x0 - x1\n        orig_dtype = self.diff.dtype\n        dtype = _get_intermediate_dtype(orig_dtype)\n        diff = self.diff.ravel().astype(dtype, copy=False)\n        return numpy.array(abs(diff).sum() / diff.size, dtype=orig_dtype),\n\n    def forward_gpu(self, inputs):\n        x0, x1 = inputs\n        self.diff = x0 - x1\n        orig_dtype = self.diff.dtype\n        dtype = _get_intermediate_dtype(orig_dtype)\n        diff = self.diff.ravel().astype(dtype, copy=False)\n        return (abs(diff).sum() / diff.dtype.type(diff.size)).astype(\n            orig_dtype, copy=False),\n\n    def backward(self, indexes, grad_outputs):\n        gy, = grad_outputs\n        coeff = gy * gy.data.dtype.type(1. / self.diff.size)\n        coeff = chainer.functions.broadcast_to(coeff, self.diff.shape)\n        gx0 = coeff * backend.get_array_module(gy.data).sign(self.diff)\n        return gx0, -gx0\n\n\ndef mean_absolute_error(x0, x1):\n    \"\"\"Mean absolute error function.\n\n    The function computes the mean absolute error between two variables. The\n    mean is taken over the minibatch. Args ``x0`` and ``x1`` must have the\n    same dimensions. This function first calculates the absolute value\n    differences between the corresponding elements in x0 and x1, and then\n    returns the mean of those differences.\n\n    Args:\n        x0 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.\n        x1 (:class:`~chainer.Variable` or :ref:`ndarray`): Input variable.\n\n    Returns:\n        ~chainer.Variable:\n            A variable holding an array representing the mean absolute\n            error of two inputs.\n\n    .. admonition:: Example\n\n        1D array examples:\n\n        >>> x = np.array([1, 2, 3]).astype(np.float32)\n        >>> y = np.array([0, 0, 0]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(2.)\n        >>> x = np.array([1, 2, 3, 4, 5, 6]).astype(np.float32)\n        >>> y = np.array([7, 8, 9, 10, 11, 12]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(6.)\n\n        2D array example:\n\n        In this example, there are 4 elements, and thus 4 errors\n        >>> x = np.array([[1, 2], [3, 4]]).astype(np.float32)\n        >>> y = np.array([[8, 8], [8, 8]]).astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(5.5)\n\n        3D array example:\n\n        In this example, there are 8 elements, and thus 8 errors\n        >>> x = np.reshape(np.array([1, 2, 3, 4, 5, 6, 7, 8]), (2, 2, 2))\n        >>> y = np.reshape(np.array([8, 8, 8, 8, 8, 8, 8, 8]), (2, 2, 2))\n        >>> x = x.astype(np.float32)\n        >>> y = y.astype(np.float32)\n        >>> F.mean_absolute_error(x, y)\n        variable(3.5)\n\n    \"\"\"\n    return MeanAbsoluteError().apply((x0, x1))[0]\n", "path": "chainer/functions/loss/mean_absolute_error.py"}]}