EE263 Homework 4

课程主页:https://see.stanford.edu/Course/EE263

这次回顾EE263作业4。

5.2

(a)

因为

所以

(b)

另解:对(1)中取

那么由于

即虚部为$0​$,所以

(c)

(d)

(e)补充证明如下结论:

首先,我们有

所以

其次我们有

所以

因此结论成立。

由之前的记号可得

由最小二乘法,我们知道该问题的解为

下面证明

证明如下:

6.2

(a)首先目标函数为

$\forall k, j-1< k \le j $,那么由牛顿第二定律可得,

任取整数$k​$,我们有

所以题目中的约束条件为

写成矩阵的形式为

其中

所以优化问题为

构造拉格朗日乘子:

求梯度并令梯度为$0$得到:

解得

matlab中使用如下命令计算即可:

1
x = pinv(A) * b

完整代码如下:

1
2
3
4
5
6
7
8
% (a)
a1 = linspace(9.5, 0.5, 10);
a2 = ones(1,10);
a3 = [linspace(4.5, 0.5, 5), 0 0 0 0 0];

A = [a1; a2; a3];
b = [1; 0; 0];
x = pinv(A) * b
1
2
3
4
5
6
7
8
9
10
11
x =
-0.0455
-0.0076
0.0303
0.0682
0.1061
0.0939
0.0318
-0.0303
-0.0924
-0.1545

作图代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
% 作图
% f(t)
figure(1);
subplot(3,1,1);
stairs(0:10, [x; 0]);
grid on;
xlabel('t');
ylabel('f(t)');
axis([0, 10, min(x) - 0.1, max(x) + 0.1]);

% p_dot(t)
T1 = toeplitz(ones(10,1), [1,zeros(1,9)]);
p_dot = T1 * x;
subplot(3,1,2);
plot(linspace(0, 10, 10), p_dot);
grid on;
xlabel('t');
ylabel('p_{dot}(t)');
axis([0, 10, min(p_dot) - 0.1, max(p_dot) + 0.1]);

% p(t)
T2 = toeplitz(linspace(0.5, 9.5, 10)', [0.5, zeros(1,9)]);
p = T2 * x;
subplot(3,1,3);
plot(linspace(0, 10, 10), p);
grid on;
xlabel('t');
ylabel('p(t)');
axis([0, 10, min(p) - 0.1, max(p) + 0.1]);

结果如下:

(b)首先

其次$\forall k, j-1< k \le j $,那么由牛顿第二定律可得,

以及任取整数$k$:

此时的约束条件为

那么

所以优化函数为

因此

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% (b)
% 作图
N = 50;
[d, m] = size(A);
Mu = logspace(-5, 2, N);
J1 = zeros(N, 1);
J2 = zeros(N, 1);

for i = 1:N
mu = Mu(i);
%x = inv(A' * A + mu * eye(m)) * A' * b;
x = (A' * A + mu * eye(m)) \ (A' * b);
j1 = norm(A * x - b) ^ 2;
j2 = norm(x) ^ 2;
J1(i) = j1;
J2(i) = j2;
end

figure(2);
plot(J1, J2);
axis tight;

结果如下:

6.5

将$f(x)$的定义代入

得到

所以线性方程组为

要使得该方程有解,应该找到最小的$m$,使得

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
% 初始化
m = 1;
X1 = x;

while 1
X1 = [X1, x .* X1(:, m)];
X2 = - diag(y) * X1;
X = [ones(N, 1), X1, X2];
m = m + 1;
% 判断
res = (rank(X) == rank([X, y]));

if res
break
end
end

m
1
5

计算系数:

1
2
3
4
%计算结果
c = X \ y;
a = c(1: m + 1)
b = c(m + 2: 2 * m + 1)
1
2
3
4
5
6
7
8
9
10
11
12
13
a =
0.2742
1.0291
1.2906
-5.8763
-2.6738
6.6845
b =
-1.2513
-6.5107
3.2754
17.3797
6.6845

6.12

所以

所以模型为

注意$\alpha ,\beta ​$未知,所以记

模型修改为

利用最小二乘法求解该问题即可,正规方程为

如果$\tilde A$满秩,那么该问题的解为

6.14

(a)

接着利用如下公式求梯度:

我们有

令梯度为$0​$得到

(b)实现上述算法:

1
A = Y * X' / (X * X')
1
2
3
4
5
6
7
8
9
10
11
A =
2.0299 5.0208 5.0104
0.0114 6.9999 1.0106
7.0424 -0.0025 6.9448
6.9977 3.9759 4.0024
9.0130 1.0449 6.9980
4.0119 3.9649 9.0267
4.9871 6.9723 8.0336
7.9425 6.0875 3.0174
0.0094 8.9722 -0.0385
1.0612 8.0208 7.0285

6.26

(a)题目的意思是利用$q_i,p_i$估计$\alpha ,d$,题目中的记号有些重复,这里假设

其中

假设

s

注意$a,d$未知,并且$d$的模长为$1$,所以可以设

那么模型为

估计$x$后,利用下式计算$a,d$即可:

(b)代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% Data for beam estimation problem
m = 5;
alpha = 0.5;
det_az =[ 3 10 80 150 275];
det_el =[ 88 34 30 20 50];
p =[ 1.58 1.50 2.47 1.10 0.001];

q1 = [cosd(det_el'), cosd(det_el'), sind(det_el')];
q2 = [cosd(det_az'), sind(det_az'), ones(m, 1)];
q = q1 .* q2;

x = (alpha * q) \ p';
a = norm(x)
d = x / a;

elevation = asind(d(3))
azimuth = asind(d(2) / cosd(elevation))
1
2
3
4
5
6
a =
5.0107
elevation =
38.7174
azimuth =
77.6623

7.3

(a)首先考虑特殊情形。

如果$\mu=0$,那么取

如果$\mu \to \infty$,那么取

所以此时的优化问题为

其中

因此

接着考虑一般情形。

对均方曲率进行处理,令

那么

所以

所以目标函数为

对比课件中的标准形式:

我们有

所以问题的解为

(b)对不同的$\mu $作图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
% 矩阵D
D1 = [eye(n - 2), zeros(n - 2, 2)];
D2 = [zeros(n - 2, 1), eye(n - 2), zeros(n - 2, 1)];
D3 = [zeros(n - 2, 2), eye(n - 2)];
D = D1 + D3 - 2 * D2;

% 作出不同的mu
% mu = 0
mu = 0;
g1 = (eye(n) + mu * n ^ 5 / (n - 2) * (D' * D)) \ f';

% mu = 1e-10
mu = 1e-10;
g2 = (eye(n) + mu * n ^ 5 / (n - 2) * (D' * D)) \ f';

% mu = 1e-3
mu = 1e-3;
g3 = (eye(n) + mu * n ^ 5 / (n - 2) * (D' * D)) \ f';

% mu = infty
A = [(1:n)', ones(n, 1)];
coef = A \ f';
g4 = A * coef;

% 作图
figure(1);
hold;
plot(f, '*');
plot(g1, '--');
plot(g2, '-.');
plot(g3, '-');
plot(g4, '--rs');
axis tight;
legend('f','u = 0', 'u = 10e-10', 'u = 10e-3', 'u = infinity', 'location', 'SouthEast');

作出权衡曲线:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% 权衡曲线
m = 30;
Mu = logspace(-8, -3, m);
x = zeros(m+2, 1);
y = zeros(m+2, 1);

x(1) = norm(f' - g1) ^ 2 / n;
y(1) = norm(D * g1) ^ 2 * n ^ 4 / (n - 2);
for i = 1: m
mu = Mu(i);
g = (eye(n) + mu * n ^ 5 / (n - 2) * (D' * D)) \ f';
x(i+1) = norm(f' - g) ^ 2 / n;
y(i+1) = norm(D * g) ^ 2 * n ^ 4 / (n - 2);
end

x(m+2) = norm(f' - g4) ^ 2 / n;
y(m+2) = norm(D * g4) ^ 2 * n ^ 4 / (n - 2);

figure(2)
plot(x, y);

8.2

题目的要求等价于求

转置后得到

对于例子中的情形,我们需要求

求解该线性方程组即可,由于方程数量大于未知数数量,所以使用左逆即可

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
G = [2 3; 1 0; 0 4; 1 1; -1 2];
G_tilde = [-3 -1; -1 0; 2 -3; -1 -3; 1 2];

A = [G'; G_tilde'];
b = [eye(2); eye(2)];

X = pinv(A) * b;
h = X';

%验证结果
h * G
h * G_tilde
1
2
3
4
5
6
ans =
1.0000 -0.0000
-0.0000 1.0000
ans =
1.0000 0
0.0000 1.0000

本文标题:EE263 Homework 4

文章作者:Doraemonzzz

发布时间:2019年06月01日 - 11:37:00

最后更新:2019年06月01日 - 11:41:30

原始链接:http://doraemonzzz.com/2019/06/01/EE263 Homework 4/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。