V2EX = way to explore
V2EX 是一个关于分享和探索的地方
现在注册
已注册用户请  登录
V2EX 提问指南
oIMOo
V2EX  ›  问与答

求助矩阵方面的计算知识

  •  
  •   oIMOo · 2019-03-11 22:45:41 +08:00 · 1506 次点击
    这是一个创建于 2090 天前的主题,其中的信息可能已经有所发展或是发生改变。

    已知以下两个值:

    • A * A^T 记为 alpha
    • A^T * A 记为 beta

    注:A 为矩阵,A^T 为 A 的转制

    求:A

    提示:SVD 奇异值分解

    实在摸不到头脑,求大佬讲解,甚至关键词都可以。

    第 1 条附言  ·  2019-03-11 23:58:49 +08:00
    import numpy as np

    A = np.array([[2, 4],[1,3],[0,0],[0,0]])
    AT = np.transpose(A)

    SDu = np.matmul(A, AT) #alpha
    SDv = np.matmul(AT, A) #beta

    UA, ZA, VA = np.linalg.svd(A)
    UAT, ZAT, VAT = np.linalg.svd(AT)

    Uu, Zu, Vu = np.linalg.svd(SDu)
    Uv, Zv, Vv = np.linalg.svd(SDv)
    第 2 条附言  ·  2019-03-12 00:05:48 +08:00
    我的代码错了么???

    np.matmul(np.matmul(Uu, Zu), Vu) != SDu

    不应该啊……
    第 3 条附言  ·  2019-03-12 00:16:37 +08:00
    https://blog.csdn.net/u012162613/article/details/42214205

    linalg.svd() 的返回值省略了一部分……
    我现在去找怎么恢复……
    第 4 条附言  ·  2019-03-12 00:22:04 +08:00
    因为 SDu 是 4*4 的矩阵,所以恢复 Sigma_u:

    Su = np.zeros([4,4])
    for i in range(4):
    Su[i][i] = Zu[i]


    Sv = np.zeros([2,2])
    for i in range(2):
    Sv[i][i] = Zu[i]

    这样确实 np.matmul(np.matmul(Uu, Su), Vu) = SDu 了
    第 5 条附言  ·  2019-03-12 09:19:39 +08:00

    这是最新的代码。

    然而注释掉的 A 可以算出正确的答案。 未注释的却不行。

    
    import numpy as np
    
    # A = np.array([[2, 4], [1, 3], [0, 0], [0, 0]])
    A = np.array([[641, 28], [256, 621], [625, 800]])
    # A = np.random.randint(0, 999, size=[3, 2])
    AT = np.transpose(A)
    
    SDu = np.dot(A, AT)  # alpha
    SDv = np.dot(AT, A)  # beta
    
    Uu, Sigmau, Vu = np.linalg.svd(SDu)
    Uv, Sigmav, Vv = np.linalg.svd(SDv)
    
    size_u,  = Sigmau.shape
    size_v,  = Sigmav.shape
    
    if(size_v > size_u):
        Sigma = np.hstack(
            (np.diag(np.sqrt(Sigmau)), np.zeros((size_v, size_v - size_u))))
    else:
        Sigma = np.vstack(
            (np.diag(np.sqrt(Sigmav)), np.zeros((size_u - size_v, size_v))))
    
    res_A = np.dot(np.dot(Uu, Sigma), Vv)
    
    print('A:\n%s\nAnswer:\n%s' % (A, res_A))
    
    9 条回复    2019-03-12 18:55:54 +08:00
    xml123
        1
    xml123  
       2019-03-11 23:06:47 +08:00
    不是有提示了吗……就是奇异值分解的过程啊,你找本矩阵论的书看一下奇异值分解的计算过程就知道了。
    oIMOo
        2
    oIMOo  
    OP
       2019-03-11 23:35:00 +08:00
    我找到了这篇文章:
    https://liam.page/2017/11/22/SVD-for-Human-Beings/
    其中 “如何计算 SVD ” *可能*与我的问题相关。

    理想情况是:
    针对 A 可以分解成 UΣV。

    根据上述理想情况:
    - A * A^T 可以分别分解成 UΣV * VΣU.
    - A^T * A 可以分别分解成 VΣU * UΣV.

    所以,还是不知道怎么求 A 啊……

    相关代码:


    import numpy as np

    A = np.array([[2, 4],[1,3],[0,0],[0,0]])
    AT = np.transpose(A)

    SDu = np.matmul(A, AT) #alpha
    SDv = np.matmul(AT, A) #beta

    UA, ZA, VA = np.linalg.svd(A)
    UAT, ZAT, VAT = np.linalg.svd(AT)

    Uu, Zu, Zu = np.linalg.svd(SDu)
    Uv, Zv, Zv = np.linalg.svd(SDv)

    -------

    另外有一篇论文。
    说到:
    - A * A^T = (U Σ V^T) * (V Σ^T U^T) = U Z Z^T U^T
    - A^T * A = (V Σ^T U^T) * (U Σ V^T) = V Z^T Z V^T

    同样使用上面的代码,可以验证这个式子。
    但是,我怎么求 A ???
    也就是说,让我只有两个乘积( alpha、beta )的前提下,我怎么把他们拆开???
    geelaw
        3
    geelaw  
       2019-03-11 23:49:03 +08:00 via iPhone
    这个问题不完整,你显然是不可能知道 A 具体是谁的——比如 A 等于正负 1 时你拿到的两个矩阵相同。

    但想要求一个 A 并不困难,你可以选 UDV',其中 U、V 分别正交相似对角化两个你知道的矩阵,D 是奇异值。
    oIMOo
        4
    oIMOo  
    OP
       2019-03-11 23:53:43 +08:00
    @geelaw 谢谢您的回复。

    您能具体说一下么?

    对 alpha 求出对应的三个值,Uu,Zu (也就是 sigma ),Vu
    对 beta 求出对应的三个值,Uv,Zv (也就是 sigma ),Vv

    然后具体怎么计算 A 呢?
    我卡在这个逻辑点了……

    谢谢
    oIMOo
        5
    oIMOo  
    OP
       2019-03-11 23:56:46 +08:00
    上面代码错了

    import numpy as np

    A = np.array([[2, 4],[1,3],[0,0],[0,0]])
    AT = np.transpose(A)

    SDu = np.matmul(A, AT) #alpha
    SDv = np.matmul(AT, A) #beta

    UA, ZA, VA = np.linalg.svd(A)
    UAT, ZAT, VAT = np.linalg.svd(AT)

    Uu, Zu, Vu = np.linalg.svd(SDu)
    Uv, Zv, Vv = np.linalg.svd(SDv)
    lovestudykid
        6
    lovestudykid  
       2019-03-12 03:09:53 +08:00
    对 alpha 和 beta 分别正交对角化,然后特征值矩阵平方,特征向量矩阵分别作为奇异值分解的左右奇异向量就完事了。对角化的过程用 svd 或者 eig 都行
    oIMOo
        7
    oIMOo  
    OP
       2019-03-12 09:24:18 +08:00
    @xml123 @geelaw @lovestudykid

    你们看最后一个 append,我真的不知道错在哪里了……
    关键是如果随机生成 A 的话,有时候对,有时候不对,疯掉了……
    geelaw
        8
    geelaw  
       2019-03-12 09:55:14 +08:00 via iPhone
    @oIMOo #7 你可以读一下下#3 的前半段。
    oIMOo
        9
    oIMOo  
    OP
       2019-03-12 18:55:54 +08:00
    @geelaw
    我在测试的时候手动排除了这点,后续再考虑。
    还有人说 A 一定要是个方阵,所以关于 A 长宽不同我也可以后续再考虑。
    然而,当 A = np.array([[5, 4, 0], [3, 7, 8], [2, 8, 3]]) 也是不能计算出结果,我就一头雾水……
    关于   ·   帮助文档   ·   博客   ·   API   ·   FAQ   ·   实用小工具   ·   2639 人在线   最高记录 6679   ·     Select Language
    创意工作者们的社区
    World is powered by solitude
    VERSION: 3.9.8.5 · 24ms · UTC 10:55 · PVG 18:55 · LAX 02:55 · JFK 05:55
    Developed with CodeLauncher
    ♥ Do have faith in what you're doing.