import numpy as np
from scipy.optimize import minimize_scalar, shgo, minimize
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution, dual_annealing
minimize_scalar
res = minimize_scalar(f, bounds=(0, 2*np.pi), method="bounded")
print(res.x)
print(res.fun)
from scipy.stats import ortho_group
def make_full_rank_matrix(d, seed=0):
"""
>>> make_full_rank_matrix(3)
array([[ 0.59537814, 0.03454805, -0.05188665],
[ 0.03454805, 0.6643958 , -0.04420196],
[-0.05188665, -0.04420196, 0.60699231]])
"""
np.random.seed(seed)
D = np.diag(np.random.rand(d))
Q = ortho_group.rvs(dim=d)
A = Q.T@D@Q
return A
# test for whether the matrix made by the above method is full rank or not
num_iter = 10
success = 0
for seed in range(num_iter):
d = 500
S = make_full_rank_matrix(d, seed)
if np.linalg.matrix_rank(S) == d:
success += 1
print(success)
別の手法
次の手法で生成される行列を考える:
$n\times k$ 次元の行列 $W$ をランダムに生成する
対角成分が正であるような対角行列 $D$ をランダムに生成する
行列 $B$ を $B = WW^\top + D$ で生成する
行列 $B$ を相似変換して行列 $A$ を得る，i.e., $A = \text{diag}[1/\sqrt{B_{ii}}] B \text{diag}[1/\sqrt{B_{ii}}]$.
実装:
def make_full_rank_matrix(d, k, seed=0):
"""
>>> make_full_rank_matrix(3, 2)
array([[1. , 0.62140315, 0.74064102],
[0.62140315, 1. , 0.70247127],
[0.74064102, 0.70247127, 1. ]])
"""
np.random.seed(seed)
W = np.random.rand(d, k)
B = W@W.T + np.diag(np.random.rand(d))
A = np.diag(1/np.sqrt(np.diag(B))) @ B @ np.diag(1/np.sqrt(np.diag(B)))
return A
# test for whether the matrix made by the above method is full rank or not
num_iter = 10
success = 0
for seed in range(num_iter):
d = 500
S = make_full_rank_matrix(d, seed)
if np.linalg.matrix_rank(S) == d:
success += 1
print(success)
x = cp.Variable(shape=n)
f = cp.norm(x, 1)
# Solve with CVXPY.
cp.Problem(cp.Minimize(f), [A@x == b]).solve()
print("Optimal value from CVXPY: {}".format(f.value))
mu = 10.0
resid = A@x - b
# penalty function method
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
print("Optimal value from method of penalty function method: {}".format(f.value))
実行結果は次のようになります．
Optimal value from method of penalty function method: 2.693372134825692
mu = 200.0 # changed
resid = A@x - b
# penalty function method
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
print("Optimal value from method of penalty function method: {}".format(f.value))
このときの実行結果は次のようになります．
Optimal value from method of penalty function method: 2.8040009637834995
f_opts = []
mus = [50.0 * i for i in range(1, 21)]
for mu in mus:
x = cp.Variable(shape=n)
f = cp.norm(x, 1)
resid = A@x - b
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
f_opts.append(f.value)
print("Optimal value from method of penalty function method: {}".format(f.value))
