diff --git a/test/test_transforms.py b/test/test_transforms.py
index ec514fe6475..1bbe1165f93 100644
--- a/test/test_transforms.py
+++ b/test/test_transforms.py
@@ -1098,16 +1098,37 @@ def _to_3x3_inv(inv_result_matrix):
         def _test_transformation(a, t, s, sh):
             a_rad = math.radians(a)
             s_rad = [math.radians(sh_) for sh_ in sh]
+            cx, cy = cnt
+            tx, ty = t
+            sx, sy = s_rad
+            rot = a_rad
+
             # 1) Check transformation matrix:
-            c_matrix = np.array([[1.0, 0.0, cnt[0]], [0.0, 1.0, cnt[1]], [0.0, 0.0, 1.0]])
-            c_inv_matrix = np.linalg.inv(c_matrix)
-            t_matrix = np.array([[1.0, 0.0, t[0]],
-                                 [0.0, 1.0, t[1]],
-                                 [0.0, 0.0, 1.0]])
-            r_matrix = np.array([[s * math.cos(a_rad + s_rad[1]), -s * math.sin(a_rad + s_rad[0]), 0.0],
-                                 [s * math.sin(a_rad + s_rad[1]), s * math.cos(a_rad + s_rad[0]), 0.0],
-                                 [0.0, 0.0, 1.0]])
-            true_matrix = np.dot(t_matrix, np.dot(c_matrix, np.dot(r_matrix, c_inv_matrix)))
+            C = np.array([[1, 0, cx],
+                          [0, 1, cy],
+                          [0, 0, 1]])
+            T = np.array([[1, 0, tx],
+                          [0, 1, ty],
+                          [0, 0, 1]])
+            Cinv = np.linalg.inv(C)
+
+            RS = np.array(
+                [[s * math.cos(rot), -s * math.sin(rot), 0],
+                 [s * math.sin(rot), s * math.cos(rot), 0],
+                 [0, 0, 1]])
+
+            SHx = np.array([[1, -math.tan(sx), 0],
+                            [0, 1, 0],
+                            [0, 0, 1]])
+
+            SHy = np.array([[1, 0, 0],
+                            [-math.tan(sy), 1, 0],
+                            [0, 0, 1]])
+
+            RSS = np.matmul(RS, np.matmul(SHy, SHx))
+
+            true_matrix = np.matmul(T, np.matmul(C, np.matmul(RSS, Cinv)))
+
             result_matrix = _to_3x3_inv(F._get_inverse_affine_matrix(center=cnt, angle=a,
                                                                      translate=t, scale=s, shear=sh))
             self.assertLess(np.sum(np.abs(true_matrix - result_matrix)), 1e-10)
diff --git a/torchvision/transforms/functional.py b/torchvision/transforms/functional.py
index a8fdbef86bf..8ae75f84c5b 100644
--- a/torchvision/transforms/functional.py
+++ b/torchvision/transforms/functional.py
@@ -8,6 +8,7 @@
 except ImportError:
     accimage = None
 import numpy as np
+from numpy import sin, cos, tan
 import numbers
 import collections
 import warnings
@@ -736,40 +737,52 @@ def _get_inverse_affine_matrix(center, angle, translate, scale, shear):
     # where T is translation matrix: [1, 0, tx | 0, 1, ty | 0, 0, 1]
     #       C is translation matrix to keep center: [1, 0, cx | 0, 1, cy | 0, 0, 1]
     #       RSS is rotation with scale and shear matrix
-    #       RSS(a, scale, shear) = [ cos(a + shear_y)*scale    -sin(a + shear_x)*scale     0]
-    #                              [ sin(a + shear_y)*scale    cos(a + shear_x)*scale     0]
-    #                              [     0                  0          1]
+    #       RSS(a, s, (sx, sy)) =
+    #       = R(a) * S(s) * SHy(sy) * SHx(sx)
+    #       = [ s*cos(a - sy)/cos(sy), s*(-cos(a - sy)*tan(x)/cos(y) - sin(a)), 0 ]
+    #         [ s*sin(a + sy)/cos(sy), s*(-sin(a - sy)*tan(x)/cos(y) + cos(a)), 0 ]
+    #         [ 0                    , 0                                      , 1 ]
+    #
+    # where R is a rotation matrix, S is a scaling matrix, and SHx and SHy are the shears:
+    # SHx(s) = [1, -tan(s)] and SHy(s) = [1      , 0]
+    #          [0, 1      ]              [-tan(s), 1]
+    #
     # Thus, the inverse is M^-1 = C * RSS^-1 * C^-1 * T^-1
 
-    angle = math.radians(angle)
-    if isinstance(shear, (tuple, list)) and len(shear) == 2:
-        shear = [math.radians(s) for s in shear]
-    elif isinstance(shear, numbers.Number):
-        shear = math.radians(shear)
+    if isinstance(shear, numbers.Number):
         shear = [shear, 0]
-    else:
+
+    if not isinstance(shear, (tuple, list)) and len(shear) == 2:
         raise ValueError(
             "Shear should be a single value or a tuple/list containing " +
             "two values. Got {}".format(shear))
-    scale = 1.0 / scale
+
+    rot = math.radians(angle)
+    sx, sy = [math.radians(s) for s in shear]
+
+    cx, cy = center
+    tx, ty = translate
+
+    # RSS without scaling
+    a = cos(rot - sy) / cos(sy)
+    b = -cos(rot - sy) * tan(sx) / cos(sy) - sin(rot)
+    c = sin(rot - sy) / cos(sy)
+    d = -sin(rot - sy) * tan(sx) / cos(sy) + cos(rot)
 
     # Inverted rotation matrix with scale and shear
-    d = math.cos(angle + shear[0]) * math.cos(angle + shear[1]) + \
-        math.sin(angle + shear[0]) * math.sin(angle + shear[1])
-    matrix = [
-        math.cos(angle + shear[0]), math.sin(angle + shear[0]), 0,
-        -math.sin(angle + shear[1]), math.cos(angle + shear[1]), 0
-    ]
-    matrix = [scale / d * m for m in matrix]
+    # det([[a, b], [c, d]]) == 1, since det(rotation) = 1 and det(shear) = 1
+    M = [d, -b, 0,
+         -c, a, 0]
+    M = [x / scale for x in M]
 
     # Apply inverse of translation and of center translation: RSS^-1 * C^-1 * T^-1
-    matrix[2] += matrix[0] * (-center[0] - translate[0]) + matrix[1] * (-center[1] - translate[1])
-    matrix[5] += matrix[3] * (-center[0] - translate[0]) + matrix[4] * (-center[1] - translate[1])
+    M[2] += M[0] * (-cx - tx) + M[1] * (-cy - ty)
+    M[5] += M[3] * (-cx - tx) + M[4] * (-cy - ty)
 
     # Apply center translation: C * RSS^-1 * C^-1 * T^-1
-    matrix[2] += center[0]
-    matrix[5] += center[1]
-    return matrix
+    M[2] += cx
+    M[5] += cy
+    return M
 
 
 def affine(img, angle, translate, scale, shear, resample=0, fillcolor=None):