2

方程式は

$C'_{ijkl} = Q_{im} Q_{jn} C_{mnop} (Q^{-1})_{ok} (Q^{-1})_{pl}$

使えました

np.einsum('im,jn,mnop,ok,pl', Q, Q, C, Q_inv, Q_inv)

仕事をすること、また期待すること

np.tensordot(np.tensordot(np.tensordot(Q, np.tensordot(Q, C, axes=[1,1]), axes=[1,0]), Q_inv, axes=[2,0]), Q_inv, axes=[3,0])

動作しますが、動作しません。

仕様:

C は 4 次弾性テンソルです。

array([[[[ 552.62389047,   -0.28689554,   -0.32194701],
         [  -0.28689554,  118.89168597,   -0.65559912],
         [  -0.32194701,   -0.65559912,  130.21758722]],

        [[  -0.28689554,  166.02923119,   -0.00000123],
         [ 166.02923119,    0.49494431,   -0.00000127],
         [  -0.00000123,   -0.00000127,   -0.57156702]],

        [[  -0.32194701,   -0.00000123,  165.99413061],
         [  -0.00000123,   -0.64666809,   -0.0000013 ],
         [ 165.99413061,   -0.0000013 ,    0.42997465]]],


       [[[  -0.28689554,  166.02923119,   -0.00000123],
         [ 166.02923119,    0.49494431,   -0.00000127],
         [  -0.00000123,   -0.00000127,   -0.57156702]],

        [[ 118.89168597,    0.49494431,   -0.64666809],
         [   0.49494431,  516.15898907,   -0.33132485],
         [  -0.64666809,   -0.33132485,  140.09010389]],

        [[  -0.65559912,   -0.00000127,   -0.0000013 ],
         [  -0.00000127,   -0.33132485,  165.98553869],
         [  -0.0000013 ,  165.98553869,    0.41913346]]],


       [[[  -0.32194701,   -0.00000123,  165.99413061],
         [  -0.00000123,   -0.64666809,   -0.0000013 ],
         [ 165.99413061,   -0.0000013 ,    0.42997465]],

        [[  -0.65559912,   -0.00000127,   -0.0000013 ],
         [  -0.00000127,   -0.33132485,  165.98553869],
         [  -0.0000013 ,  165.98553869,    0.41913346]],

        [[ 130.21758722,   -0.57156702,    0.42997465],
         [  -0.57156702,  140.09010389,    0.41913346],
         [   0.42997465,    0.41913346,  486.62412063]]]])

Q は、x 座標と y 座標を変更する回転行列です。

array([[ 0,  1,  0],
       [-1,  0,  0],
       [ 0,  0,  1]])

Q_inv は

array([[-0., -1., -0.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.]])

np.einsumにつながる

array([[[[ 516.15898907,   -0.49494431,   -0.33132485],
         [  -0.49494431,  118.89168597,    0.64666809],
         [  -0.33132485,    0.64666809,  140.09010389]],

        [[  -0.49494431,  166.02923119,    0.00000127],
         [ 166.02923119,    0.28689554,   -0.00000123],
         [   0.00000127,   -0.00000123,    0.57156702]],

        [[  -0.33132485,    0.00000127,  165.98553869],
         [   0.00000127,   -0.65559912,    0.0000013 ],
         [ 165.98553869,    0.0000013 ,    0.41913346]]],


       [[[  -0.49494431,  166.02923119,    0.00000127],
         [ 166.02923119,    0.28689554,   -0.00000123],
         [   0.00000127,   -0.00000123,    0.57156702]],

        [[ 118.89168597,    0.28689554,   -0.65559912],
         [   0.28689554,  552.62389047,    0.32194701],
         [  -0.65559912,    0.32194701,  130.21758722]],

        [[   0.64666809,   -0.00000123,    0.0000013 ],
         [  -0.00000123,    0.32194701,  165.99413061],
         [   0.0000013 ,  165.99413061,   -0.42997465]]],


       [[[  -0.33132485,    0.00000127,  165.98553869],
         [   0.00000127,   -0.65559912,    0.0000013 ],
         [ 165.98553869,    0.0000013 ,    0.41913346]],

        [[   0.64666809,   -0.00000123,    0.0000013 ],
         [  -0.00000123,    0.32194701,  165.99413061],
         [   0.0000013 ,  165.99413061,   -0.42997465]],

        [[ 140.09010389,    0.57156702,    0.41913346],
         [   0.57156702,  130.21758722,   -0.42997465],
         [   0.41913346,   -0.42997465,  486.62412063]]]])

これは正しいと思いますが、4 つnp.tensordotは

array([[[[ 552.62389047,   -0.28689554,    0.32194701],
         [  -0.28689554,  118.89168597,    0.65559912],
         [  -0.32194701,   -0.65559912, -130.21758722]],

        [[  -0.28689554,  166.02923119,    0.00000123],
         [ 166.02923119,    0.49494431,    0.00000127],
         [  -0.00000123,   -0.00000127,    0.57156702]],

        [[  -0.32194701,   -0.00000123, -165.99413061],
         [  -0.00000123,   -0.64666809,    0.0000013 ],
         [ 165.99413061,   -0.0000013 ,   -0.42997465]]],


       [[[  -0.28689554,  166.02923119,    0.00000123],
         [ 166.02923119,    0.49494431,    0.00000127],
         [  -0.00000123,   -0.00000127,    0.57156702]],

        [[ 118.89168597,    0.49494431,    0.64666809],
         [   0.49494431,  516.15898907,    0.33132485],
         [  -0.64666809,   -0.33132485, -140.09010389]],

        [[  -0.65559912,   -0.00000127,    0.0000013 ],
         [  -0.00000127,   -0.33132485, -165.98553869],
         [  -0.0000013 ,  165.98553869,   -0.41913346]]],


       [[[   0.32194701,    0.00000123,  165.99413061],
         [   0.00000123,    0.64666809,   -0.0000013 ],
         [-165.99413061,    0.0000013 ,    0.42997465]],

        [[   0.65559912,    0.00000127,   -0.0000013 ],
         [   0.00000127,    0.33132485,  165.98553869],
         [   0.0000013 , -165.98553869,    0.41913346]],

        [[-130.21758722,    0.57156702,    0.42997465],
         [   0.57156702, -140.09010389,    0.41913346],
         [  -0.42997465,   -0.41913346,  486.62412063]]]])

負の大きな数に注意してください。

4

1 に答える 1

3

アプローチ #1

1 つの方法は、単一のステップではなく、信頼できるnp.tensordot人からの助けを借りて、 と同じ結果を得るために使用することです-np.einsum broadcasting

# Get broadcasted elementwise multiplication between two versions of Q.
# This corresponds to "np.einsum('im,jn,..', Q, Q)" producing "'ijmn"" 
# broadcasted version of elementwise multiplications between Q's.
Q_ext = Q[:,None,:,None]*Q[:,None,:]

# Similarly for Q_inv : For "np.einsum('..ok,pl', Q_inv, Q_inv)" get "'opkl'" 
# broadcasted version of elementwise multiplications between Q_inv's.
Q_inv_ext = Q_inv[:,None,:,None]*Q_inv[:,None,:]

# Perform "np.einsum('im,jn,mnop,ok,pl', Q, Q, C)" with "np.tensordot".
# Notice that we are using the last two axes from 'Q_ext', so "axes=[2,3]"
# and first two from 'C', so "axes=[0,1]" for it. 
# These axes would be reduced by the dot-product, leaving us with 'ijop'.
parte1 = np.tensordot(Q_ext,C,axes=([2,3],[0,1]))

# Do it one more time to perform "np.einsum('ijop,ok,pl', parte1,Q_inv,Q_inv)"
# to reduce dimensions represented by 'o,p', leaving us with 'ijkl'.
# To confirm, compare the following against original einsum approach :
# "np.einsum('im,jn,mnop,ok,pl->ijkl', Q, Q, C, Q_inv, Q_inv)"
out = np.tensordot(parte1,Q_inv_ext,axes=([2,3],[0,1]))

アプローチ #2

broadcastingのインスタンスを 2 つ以上使用することを避けたい場合はnp.tensordot、次のようにします。

# Perform "np.einsum('jn,mnop', Q, C). Notice how, Q is represented by 'jn' 
# and C by 'mnop'. We need to reduce the 'm' dimension, i.e. reduce 'axes=1'
# from Q and `axes=1` from C corresponding to `n' in each of the inputs.
# Thus, 'jn' + 'mnop' => 'jmop' after 'n' is reduced and order is maintained.
Q_C1 = np.tensordot(Q,C,axes=([1],[1]))

# Perform "np.einsum('im,jn,mnop', Q, Q, C). We need to use Q and Q_C1.
# Q is 'im' and Q_C1 is 'jmop'. Thus, again we need to reduce 'axes=1'
# from Q and `axes=1` from Q_C1 corresponding to `m' in each of the inputs.
# Thus, 'im' + 'jmop' => 'ijop' after 'm' is reduced and order is maintained.
parte1 = np.tensordot(Q,Q_C1,axes=([1],[1]))

# Use the same philosophy to get the rest of the einsum equivalent, 
# but use parte1 and go right and use Q_inv 
out = np.tensordot(np.tensordot(parte1,Q_inv,axes=([2],[0])),Q_inv,axes=([2],[0]))

の秘訣は、パラメーターによって削減されるディメンションと、縮小されたディメンションが残りの入力のディメンションに対してどのように整列するかnp.tensordotを追跡することです。axes

于 2016-03-03T17:01:58.473 に答える