条件下における式を解きたい

Question

ある条件を満たすリストと定数を求めたい

s = [s1, s2, s3]
t = [t1, t2, t3]

が与えられた時、 以下の条件を満たす p = [p1, p2, p3] と h を求めたいです。

(条件)
・ pk = max{ tk(sk - h), 0 } k = 1,2,3
・ h は p1 + p2 + p3 = 1 を満たす為の定数

この条件下で出来るだけ p を 0 にしないようなコードを書きたいと考えております。
(0 を許可するコードなら既に書けております)

どういったコードが考えられますでしょうか。
考え方でもご教授頂けると幸いです。
よろしくお願い申し上げます。

返信が遅れまして大変申し訳ございません。
またご回答くださり誠に有難うございます。
以下、コードです。

def calculate_send_power(bits, A, r0, a, b):
    num_stream = len(bits)
    over = 0
    under = 0
    for k in range(num_stream):
        rk = A[k] * r0
        ak = a[bits[k]]
        bk = b[bits[k]]
        mk = bits[k]

        over += bk/rk * math.log(ak * mk * rk / bk)
        under += bk/rk

    xi = (over - 1) / under
    p = []

    for k in range(num_stream):
        rk = A[k] * r0
        ak = a[bits[k]]
        bk = b[bits[k]]
        mk = bits[k]

        if bk/rk * (math.log(ak * mk * rk / bk) - xi) > 0:
            p.append(bk/rk * (math.log(ak * mk * rk / bk) - xi))

    if len(p) != num_stream:
        p = calculate_send_power(bits[:-1], A, r0, a, b)

    return p

if __name__ == '__main__':
    for bits in bit_pattern:
        num_stream = len(bits)
        p = calculate_send_power(bits, A, r0, a, b)

        if len(p) != num_stream:
            for i in range(num_stream - len(p)):
                p.append(0)

かなり分かりにくいかと思いますが、 calculate_send_power() における over が s を表し、under が t を表し、 xi が h を表しているとお考えください。

以下にそれぞれの変数を述べておきます。

a = [0,0,1/2,0,3/8,0,7/24,0,15/64]
b = [0,0,2,0,10,0,42,0,170]

bit_pattern = [[8], [6, 2], [4, 4], [4, 2, 2], [2, 2, 2, 2]]

H = (4,4)の複素行列 実部、虚部ともに正規分布に従う乱数
gram_H = np.conjugate(H.T) @ H
U, A, Uh = svd(gram_H)

SNR = 0
Ps = 1.0
sigma = math.sqrt( Ps / (2 * math.pow(10.0, SNR/10)))

Pn = 2*sigma*sigma
r0 = Ps/Pn

このコードは bit_pattern = [[8], [6, 2], [4, 4], [4, 2, 2], [2, 2, 2, 2]] のそれぞれに割り当てる電力(電力総和 = 1)を決定するもので以下のようになります。

bits	p(ower)	sum
[8]	[1.000]	1.0
[6, 2]	[0.453, 0.546]	1.0

以下に実行結果を述べます。

for bits in bit_pattern:
    num_stream = len(bits)
    tmp = 0
    p = calculate_send_power(bits, A, r0, a, b)
    if len(p) != num_stream:
        for i in range(num_stream - len(p)):
            p.append(0)

    print("bits = {} , power = {}, sum = {}".format(bits, p, sum(p)))

>>> 
bits = [8] , power = [1.00], sum = 1.0
bits = [6, 2] , power = [0.39036578941718647, 0.6096342105828132], sum = 1.0
bits = [4, 4] , power = [0.6905457297074978, 0.30945427029250216], sum = 1.0
bits = [4, 2, 2] , power = [0.4747411414947003, 0.35540013052184316, 0.1698587279834565], sum = 1.0
bits = [2, 2, 2, 2] , power = [0.3157712110211322, 0.40062777763337526, 0.2836010113454926, 0], sum = 1.0

と出来ているのですが、 bits = [2, 2, 2, 2] の4項目に割り当てる電力が 0 になってしまいます。この条件下で出来るだけ p を 0 にしないようなコードを書きたいと考えており、考え方でもご教授頂けると幸いです。よろしくお願い申し上げます。

Comments

[can110]

手計算で求めることはできますか？また、その手順を説明することはできますか？

[can110]

「p に負があるなら～」や「出来る限り p には 0 を入れたくない」などの制約条件があれば質問に記載ください。
また、重ねて確認ですが
手計算で求めることはできますか？また、その手順を説明することはできますか？
あるいは現状結果を求めることができるコードは提示できますか？

[Zuishin]

解が複数ある可能性がありますね。

[sfdust]

回答のコード修正したんですが、３つともゼロでない場合を求める強い制約となっています。
（３つのpのうちいずれか1つ以上ゼロがありうる場合は即座にNoneを返す）
ただ、質問の元のコードでは、p1,p2,p3の３つではなくて、4つでやってますよね？
けれど４つともにゼロでない場合はないような気もします。
結局、求めたいのは、「４つのうちゼロであるものが最小となるすべてのケース」ということでしょうか。そうなるとアプロ―チが異なってくるのですが・・・

[wandora]

ご返信頂き誠に有難うございます。
またコードも頂き重ねて御礼申し上げます。

上のコードは仰ってくださった線形アプローチで h を解いていまして、もし p に 0 が入った場合は p[:-1] を再帰的にかけて、足りない部分に 0 を入れるといったことをしております。

分かりやすいように p の個数は 3 で質問させていただきましたが、4 以上をとることもあります。こういった場合、線形で解くと p に 0 が入ってしまうことがあるので、線形で解く以外のアプローチが無いかと解答の考え方を質問させて頂いていております。

仰ってくださった「p の個数が4つのうちゼロであるものが最小となるすべてのケース」を求める場合はどういったアプローチになりますでしょうか。

お手数お掛けして申し訳御座いません。
何卒、宜しくお願い致します。

Answer 1

s = [s[0], s[1], s[2]]
t = [t[0], t[1], t[2]]

としたとき、max()を無視すれば、

p1 = t[0]*s[0]-t[0]*h
p2 = t[1]*s[1]-t[1]*h
p3 = t[2]*s[2]-t[2]*h

となります。
p1～p3のいずれかが0となるようなhは除外するのですから、結局、
直交座標面(h, y)において
直線y = p1+p2+p3 = t[0]*s[0]-t[0]*h + t[1]*s[1]-t[1]*h + t[2]*s[2]-t[2]*h
= -(t[0] + t[1] + t[2])*h + (t[0]*s[0] + t[1]*s[1] + t[2]*s[2])

が直線 y=1と交わるときのh座標を算出し
求めたhがいずれのk=0,1,2についてもpk=tk(sk - h)>0を満たしているか
確認すればよいのではないでしょうか。

上記で
x = t[0] + t[1] + t[2]
k = t[0]*s[0] + t[1]*s[1] + t[2]*s[2]
とおけば
y = -xh+k　と　y = 1の交点を求めるわけですから、
h = (k-1)/x
となります。（x=0の場合は除外）

※s,tの中身によっては、解を持ちません。

def func1(s, t):
    x = sum(t)  # t[0] + t[1] + t[2]...
    k = sum(ss*tt for ss, tt in zip(s, t))  # [0]*s[0] + t[1]*s[1] + t[2]*s[2]...

    if x == 0:
        h = None
    else:
        h = (k-1)/x
        hh = h
        for ss, tt in zip(s, t):
            if tt*(ss - hh) <= 0:
                h = None
    if h is None:
        return None
    return tuple(tt*(ss - h) for ss, tt in zip(s, t))

実行例

s = [1, 2, 3]
t = [0.4, 0.3, 0.2]

r = func1(s,t)
print("answer=", r)
print("sum = {}".format(sum(r) if r else None))

# answer= (0.13333333333333328, 0.3999999999999999, 0.4666666666666666)
# sum = 0.9999999999999998

s = [1, 2, 3, 2]
t = [0.2, 0.3, 0.2, 0.1]
# answer= (0.04999999999999996, 0.37499999999999994, 0.45, 0.12499999999999999)
# sum = 1.0

s = [1, 2, 3, 2]
t = [0.2, 0.3, 0.2, 0.5]
# answer= None
# sum = None