叶藏也不行 | 博客

~~~


  • Home

  • About

  • Archives

  • Tags

  • Search

进程管理与调度

Posted on 2023-07-31

进程管理和调度

1 进程优先级

进程类型

  • 硬实时进程,Linux主流不支持

  • 软实时进程

  • 普通进程,可分配优先级

抢占式多任务处理;CFS,调度实体,

2 进程生命周期

进程状态

  • 运行

  • 等待

  • 睡眠

进程表,僵尸进程(资源已释放,但进程表项还存在);

抢占式多任务处理:

  • 普通进程总是可能被抢占

  • 内核态系统调用无法被其他进程抢占,但中断可以中止系统调用

  • 中断可以暂停用户态和内核态进程,具有最高优先级

内核2.5后内核抢占支持内核系统调用被其他进程抢占;

3 进程表示

struct task_struct结构

  • 状态和执行信息

  • 虚拟内存信息

  • 进程身份凭据

  • 文件和文件系统信息

  • 进程特定于CPU的运行数据线程信息

  • 进程间通信有关信息

  • 信号处理程序

state进程当前状态

  • TASK_RUNNING

  • TASK_INTERRUPTIBLE

  • TASK_UNINTERRUPTIBLE

  • TASK_STOPPED

  • TASK_TRACED

  • EXIT_ZOMBIE, EXIT_DEAD

资源限制rlimit

3.1 进程类型

  • fork

  • exec

  • clone,用于实现线程

3.2 命名空间

  1. 概念:虚拟化,容器,建立系统不同的视图,chroot机制,

    • fork或clone时可选择命名空间选项

    • unshare系统调用

  2. 实现:uts、ipc、mnt、pid、user、net等命名空间,struct nsproxy;

  3. 进程ID号:全局ID、局部ID;pid,upid,pid_namespace等之间的关系:

    • upid: (number, pid_namespace)

    • pid_namespace: number->pid

    • pid: (level, upid, …) pid->task_struct->pid

    • 会为pid的每个命名空间分配一个局部PID,使用idr(radix_tree)数据结构;

  4. 进程关系:使用链表和指针表示;

4 进程管理相关系统调用

4.1 进程复制

  • fork,COW写时复制

  • vfork,父子进程共享数据

  • clone,产生线程,对父子进程之间的共享、复制进行精确控制

kernel_clone->copy_process

4.2 内核线程

4.3 启动新程序

execve系统调用,二进制格式,

4.4 退出进程

exit系统调用,

5 调度器的实现

5.1 概观

CFS;进程等待时间红黑树排序;虚拟时钟;进程优先级Nice值;

  • 进程有不同优先级需要考虑

  • 进程不能切换太频繁,多媒体系统延迟;

5.2 数据结构

  • 调度类(实现调度策略,策略模式)

  • 任务切换(CPU相关)

task_struct的相关成员

sched_class调度器类(提供一堆回调函数)

就绪队列,各个活动进程只出现在一个就绪队列中:

  • struct rq,嵌入各种子就绪队列(cfs_rq,rt_rq,dl_rq)
  • 每个CPU包含一个就绪队列,per-cpu数据

调度实体,sched_entity类,进程是可调度实体,但可调度实体更为广义;

5.3 处理优先级

  1. 优先级的内核表示:实时优先级、普通优先级(Nice值)

  2. 计算优先级:综合考虑动态优先级、普通优先级、静态优先级;

  3. 计算负荷权重:struct load_weight,优先级权重值转换表,转换考虑实时进程(是普通进程的两倍);就绪队列也关联一个负荷权重;

5.4 核心调度器

  1. 周期性调度器scheduler_tick,按照一定周期调用该函数;

  2. 主调度器schedule:主动让出CPU给另一个进程,内核系统调用返回后检查重调度标志TIF_NEED_RESCHED;。。。

  3. 与fork的交互;

  4. 上下文切换;context_switch,switch_mm,switch_to,惰性TLB,switch_to三个参数传递两个变量prev = switch_to(prev, next), A->B->C->A, prev = C;惰性FPU模式;

6 完全公平调度类

struct fair_sched_class

6.1 数据结构

6.2 CFS操作

6.3 队列操作

6.4 选择下一个进程

6.5 处理周期性调度器

6.6 唤醒抢占

6.7 处理新进程

7 实时调度类

7.1 性质

7.2 数据结构

7.3 调度器操作

8 调度器增强

8.1 SMP调度

8.2 调度域和控制组

8.3 内核抢占和低延迟相关工作

9 小节

Read more »

进程管理与调度

Posted on 2022-09-29

《深入Linux内核架构》进程管理与调度相关笔记

进程管理相关的系统调用

  • 命令空间层次化机制,pid类型创建与管理等

  • fork, vfork, clone系统调用
  • 写时复制COW技术,父子进程共享物理内存页,标记页表,要写一页就复制一页
  • NPTL(Native POSIX Threads Library),用户空间线程库,尽管内核提供线程机制,但是用户空间需要线程库才能使用
  • 大部分工作由体系无关函数do_fork实现
  • thread_info和task_struct的关系,前者一般占两个内存页,包含内核栈以及体系结构特定的线程相关信息
  • sys_futex(),快速的用户空间互斥量;
  • 内核线程,内核函数以进程的方式并行运行;kernel_thread函数;内核线程不可修改用户空间部分内容(随机的用户空间进程内容);为强调用户空间部分不能问,mm设置为空指针。但由于内核必须知道用户空间当前包含了什么,所以在active_mm中保存了指向mm_struct的一个指针来描述它
  • execve系统调用启动新程序;可执行程序格式,ELF格式;
  • exit系统调用终止进程;编译器在main函数末尾自动添加该调用;资源释放;

调度器的实现

  • 调度策略,上下文切换;
  • 红黑树可运行队列;虚拟时间与实时时间;
  • 调度器类sched_class实现不同的调度策略;就绪队列run queue,每个CPU都有一个,一个进程只在一个rq上;
  • 调度实体sched_entity;
  • 实时互斥量问题RT-Mutex;
  • 惰性TLB;
  • 上下文切换的复杂,需要三个参数传递两个变量,因为上下文切换涉及三个进程;
  • 惰性FPU模式,内核不使用的寄存器如浮点寄存器,不会保存和恢复除非有使用;

完全公平调度类

  • 虚拟时钟的模拟;越重要的进程会有越高的优先级,更大的权重,累加的虚拟运行时间更小;
  • 在运行队列红黑树中,虚拟时钟总是增加,活动进程总是向右移动,越重要的移动越慢,因此总是倾向于保持在队列前面;睡眠进程虚拟时钟保持不变,而队列最小虚拟时钟总是单调递增,睡眠进程醒来后在红黑树中的位置更靠左(键值计算方式为虚拟时钟减去队列最小虚拟时钟),被调度的机会也变大;
  • 延迟跟踪,调度延迟,即保证每个可运行进程都至少运行一次的某个时间间隔;
  • 抢占进程,确保被强占进程执行一定的时间,不至于切换频繁;

实时调度类

  • 循环进程、先进先出进程;

调度器增强

  • SMP调度,负载均衡,亲和性,进程迁移,
  • 调度域和控制组(cgoups);
  • 内核抢占和低延迟相关工作;不仅用户空间应用程序可被中断,内核也会被中断;thread_info中的抢占计数器来防止两个处理器在同一个临界区工作;内核代码检查长时间运行的函数,并在适当之处插入cond_resched调用;
  • TASK_KILLABLE状态,进程睡眠,但不响应非致命信号,可以被致命信号杀死;
Read more »

索引优先队列

Posted on 2020-02-16

《算法 第四版》索引优先队列有感…

重点理解示意图!

一、优先队列

使用二插堆, 利用swim和sink操作很容易实现最小优先队列和最大优先队列(取决于比较方式), 同时也可以顺便实现堆排序算法, 具体细节在此按下不表, 请自行搜索学习。

二、索引优先队列

1、前言

《算法 第四版》提到了索引优先队列这个数据结构, 该数据结构不仅有优先队列的快速操作, 同时还能方便地利用索引自由地进行更多的操作, 比如返回最大(最小)元素的索引, 修改某个索引对应的元素, 删除某个索引的元素或是最大(最小)元素等等. 利用该数据结构我们可以方便地解决很多问题. 但是, 该数据结构的实现往往很容易让人迷惑, 所以, 我们用最大索引优先队列为例, 尽可能地让读者对此有清晰的理解。

2、数据结构

索引优先队列主要由三个数组实现, 其中两个索引数组, 一个元素数组(容器), 我们把它们分别叫做qp[], pq[], items[]数组. pq[]数组索引items[], qp[]数组索引pq[], 我们维护的优先队列实际上是对pq[]数组的维护. 我们首先用插入操作举个例子,先来看看它的代码:

    public void insert(int k, Item v) {
        N++;
        qp[k] = N;
        pq[N] = k;
        items[k] = v;
        swim(N);
    }

一时看不懂没关系,继续看下去,待会再返回来看代码.我们可以看到,我们对Item元素的索引k实际上就是qp[]和items[]数组的下标,对没有索引的元素,我们会将qp[]数组对应下标的值置-1表示没有这个索引值.

    public IndexMaxPQ(int max) {
        pq = new int[max + 1];
        qp = new int[max + 1];
        items = (Item[]) new Comparable[max + 1];
        for(int i = 0; i <= max; i++) {
            qp[i] = -1; // 表示该索引没有item
        }
    }

我们一旦使用了索引k,那么要么改变该索引的元素,要么删除该索引和对应的元素,否则他们是固定不变的(实际上,items[]数组不变,qp[]数组是会改变的,因为它是用来索引pq[]数组的,这点我们下面会更清楚地讲解明白). 然后我们可以看到我们使用了一个优先队列中很重要的函数swim保持其堆有序,在此处的作用也和优先队列中的一样,不过他所使用的exchange函数有很大的不同.我们再来看看下面的图,就能更清楚这三个数组之间的关系了:

在这里我们可以看到, 对于任意的1 <= k <= N,pq[qp[k]]=qp[pq[k]] = k(这也是书上的结论),qp[]和items[]数组关于pq[]数组对称. 我们使用swim和sink函数维护堆有序性的时候,实际上也是利用这种关系来实现的.维护pq[]数组时,我们利用pq[]指向的Item元素进行比较,决定是否上浮或下沉,然后重新修改pq[]和qp[]指向的元素,swim和sink函数的操作对象是pq[]数组.我们来看一看exchange函数和进行比较的less函数:

    private boolean less(int i, int j) {
        return items[pq[i]].compareTo(items[pq[j]]) < 0;
    }

    private void exch(int i, int j) {
        int k1 = pq[i];
        int k2 = pq[j];
        qp[k1] = j;
        pq[j] = k1;
        qp[k2] = i;
        pq[i] = k2;
    }

less函数就不多解释了,我们来看一看exchange函数.我们首先保存pq[]数组中i,j指向的元素的索引,因为他们是不改变的,然后修改我们qp[]数组和pq[]数组相应的索引值实现了交换操作,这个操作通过上面的示意图可以更直观的看出来. 我们一旦理解了这个过程,其他的东西便可迎刃而解了: 1、最大(最小)值或者索引,可以直接通过pq[1]得出; 2、修改某个索引元素,只需要修改items[]数组对应索引下标的元素,然后重新调整对应的pq[]数组使其堆有序(修改后的元素可能打破这个规则); 3、删除最大(最小)或者某个索引值对应元素,我们只需要像优先队列一样,将该索引指向的pq[]数组元素与最后一个元素交换后减小堆的大小,然后利用swim和sink函数调整交换后的堆使其堆有序.

3、完整代码(以最大优先队列为例)

import java.lang.Comparable;

public class IndexMaxPQ<Item extends Comparable<Item>> {
    private int[] pq;
    private int[] qp;
    private Item[] items;
    private int N = 0;

    private void swim(int k) {
        while(k > 1 && less(k/2, k)) {
            exch(k, k/2);
            k = k/2;
        }
    }

    private void sink(int k) {
        while(k*2 <= N) { // 有子节点才继续
            int j = k*2;
            if(j+1 <= N && less(j, j+1)) j += 1;

            if(!less(k, j)) break;
            exch(k, j);
            k = j;
        }
    }

    private boolean less(int i, int j) {
        return items[pq[i]].compareTo(items[pq[j]]) < 0;
    }

    private void exch(int i, int j) {
        int k1 = pq[i];
        int k2 = pq[j];
        qp[k1] = j;
        pq[j] = k1;
        qp[k2] = i;
        pq[i] = k2;
    }

    public IndexMaxPQ(int max) {
        pq = new int[max + 1];
        qp = new int[max + 1];
        items = (Item[]) new Comparable[max + 1];
        for(int i = 0; i <= max; i++) {
            qp[i] = -1; // 表示该索引没有item
        }
    }

    public boolean contains(int k) {
        return qp[k] != -1;
    }

    public void insert(int k, Item v) {
        N++;
        qp[k] = N;
        pq[N] = k;
        items[k] = v;
        swim(N);
    }

    public void change(int k, Item v) {
        items[k] = v;
        swim(qp[k]);
        sink(qp[k]);
    }

    public int delmax() {
        int k = pq[1];
        exch(1, N--);
        sink(1);
        items[pq[N+1]] = null;
        qp[pq[N+1]] = -1;
        return k;
    }

    public void delete(int k) {
        exch(qp[k], N--);
        swim(qp[k]);
        sink(qp[k]);
        items[pq[N+1]] = null;
        qp[pq[N+1]] = -1;
    }

    public int maxIndex() {
        return pq[1];
    }

    public Item max() {
        return items[pq[1]];
    }

    public boolean isEmpty() {
        return N == 0;
    }

    public int size() {
        return N;
    }

    public static void main(String[] args) {
    }
}

Read more »

K-SVD实现人脸图像恢复

Posted on 2019-12-20

本博客为表达学习课程期末实验,使用K-SVD算法实现缺失人脸图像的恢复。图像数据取自Yale B 数据集中的38张192*168的人脸正面图像…

一、实验目的(问题描述)

本实验拟使用K-SVD算法实现缺失人脸图像的恢复。图像数据取自Yale B 数据集中的38张192*168的人脸正面图像。将图像数据划分为训练集和测试集,在训练集上使用K-SVD算法得到字典,对测试集中的图像进行50%、70%的随机像素值缺失处理。使用训练集训练得到的字典填补缺失图像并计算重建人脸图像与原始图像的平均误差。

二、算法描述和分析

1、OMP(Orthogonal Matching Pursuit 正交匹配追赶法)算法

1)原理

在匹配算法(MP)基础上,每次迭代中通过在已选的集合中的全部元素上重新调整表达系数,使得每次都能最大化的减小残差。

2)性质

A、
\[rk\perp span\{a_{\lambda 1},a_{\lambda 2},a_{\lambda 3},\dots,a_{\lambda k} \}\]

#####

B、 对于m维向量,该算法最多在m步内即可收敛

3)算法步骤

2、K-SVD算法

1)字典学习

2)K-SVD原理

K-SVD是一种字典学习算法(其他算法有MOD方向法等),该算法通过两个主要步骤训练字典:

1、稀疏编码:固定字典D,求解稀疏表达矩阵X;

2、字典更新(主要创新点):在当前表达X下,逐一更新字典D中原子(列)(固定其他列不变),同时也更新表达矩阵X的对应行。

其中,第一步可以通过常见的稀疏编码算法实现(OMP,BP),本实验使用OMP算法;

第二步通过优化目标函数来更新原子(列):

由于把第k个原子剥离后,表达中产生了空洞,我们的目标就是找到一个新的dk^,更好的填补这个空洞;或者说,寻找Ek的秩为1的最优逼近,但直接求解不易,而这又可以通过核范数(University of Washington 的教授M. Fazel提出基于核范数对矩阵的秩进行封装)来进行有效估计,故可以将Ek进行SVD分解即可满足需求。

同时为了保证新生成的X的第k行gk_T还是稀疏的,我们需要对Ek和gk_T进行“压缩”,只取gkT非0元素对应下标的原子(列)和元素。

3)算法步骤

三、算法实现

1、OMP算法

1)OMP算法实现

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#-- OMP算法辅助函数
def find_max_atom_vector(rk_1, A):
    '''
        #-- 从字典A中找出最大原子向量
        #-- 输入参数:
        #-- rk_1: 前一个残差向量
        #-- A:字典
        #-- 输出:返回最大原子向量下标
    '''
    return np.argmax([np.abs(np.dot(ai,rk_1))  for i,ai in enumerate(A.T[:,]) ])


# -- OMP算法求解 b = Ax 稀疏度最小的x
def OMP(A, b, e, k):
    '''
        # -- OMP算法求解 b = Ax 稀疏度最小的x
        # -- 输入参数: 
        # -- A:  字典
        # -- b:  需要稀疏表达的向量
        # -- e:  停止条件1,当误差小于e时可以直接停止迭代
        # -- k:  停止条件2,最大迭代次数
        # -- 输出:b在A上的稀疏表达x
    '''
    N = A.shape[1]
    x = np.zeros([N])  #初始化稀疏表达
    S = []  # 支撑集
    rk = b  # 初始残差
    indexs = []  #选择的基下标

    #  迭代
    for i in range(k):
        # 选择残差投影最大的原子向量
        index = find_max_atom_vector(rk, A,)

        # 选择的下标,即最终的x不为0的位置
        indexs.append(index)

        # 加入新选择的基
        S.append(A[:, index])

        # 利用矩阵的逆重新计算表达式
        Ak = np.array(S).T
        xk = np.linalg.inv(Ak.T.dot(Ak)).dot(Ak.T.dot(b))
      
        # 更新残差
        rk = b - Ak.dot(xk)

        # 残差已经足够小了
        if(np.sqrt(np.mean(np.square(rk))) - e < 0):
            break

    # 构造最终求得的x并返回
    x[indexs] = xk.squeeze()
    return x

2)OMP算法效果展示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
m = 20
n = 10000
b = np.random.normal(size = (m, 1))
A = np.random.normal(size = (m, n))

x = OMP(A,b.copy(),1e-5,7)
plt.plot(b, 'g', label='原始数据')
plt.plot(A.dot(x), 'r', label='重构数据')

plt.xlabel('indexs')
plt.ylabel('values')
plt.title('OMP效果图')
plt.legend()

plt.show()

##### 拟合效果

结论

​ 从图中可以看出,OMP算法重构的数据对原始数据进行了较好的拟合(稀疏度为7);

3)OMP算法进一步分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
m = 20
n = 10000
b = np.random.normal(size=(m, 1))

A = np.random.normal(size=(m, n))
errors = []

for i in range(20):
    x = OMP(A, b.copy(), 1e-5, i+1)
    new_x = A.dot(x)
    error = np.sqrt(np.sum(np.square(b.squeeze()-new_x)))
    errors.append(error)

plt.plot([i+1 for i in range(20)], errors, 'g', label='误差')

plt.xlim((0,25))
plt.xlabel('稀疏度')
plt.ylabel('误差值')
plt.title('OMP效果图')
plt.legend()   
plt.show()

结论

​ 使用不同的稀疏度对原始数据进行稀疏表达后再拟合,从图中可以看出,稀疏度越大,误差越小(稀疏度从1到20)。但是,稀疏度太大的话,需要的迭代次数也更多,因此实际使用时,要从效率和效果两方面综合考虑来选择合适的稀疏度。

2、K-SVD算法

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#-- K-SVD 求解
#-- 矩阵的F范式---
def normalization_F(A):
    '''
        计算矩阵的F范式,定义见Lecture6
        输入参数:矩阵A
        输出:矩阵A的F范式
    '''
    return np.sqrt(np.sum(np.square(A)))


def K_SVD(D, S, Y, e, L ):
    '''
        利用K-SVD算法训练字典D

        输入参数:
            原始字典D:m x K 固定大小
            训练数据Y:m x N
            稀疏度S:OMP算法退出条件
            残差e:OMP算法退出条件
            迭代次数L:K-SVD迭代次数
            目标:yi = Dxi + ei
                || xi ||0 <= S <= m; || ei ||(2,2) <= e;(任意1 <= i <= N) 
        输出:
            D:训练好的字典
            X:数据Y在字典D上的系稀疏表达,大小为 K x N
    '''
    m = D.shape[0]
    K = D.shape[1]
    N = Y.shape[1]
    print('Y = DX, 其中 Y = %d x %d ; D = %d x %d ; X = %d x %d'%(m,N,m,K,K,N))

    X = np.zeros([K, N])
    for i in range(L):
        # 步骤一: 利用OMP寻找最相关稀疏矩阵X
        for i in range(N):
            X[:,i] = OMP(D,Y[:,i],e,S)

        # 步骤二:字典更新
        # 每次更新一列,固定其他列不变
        for k in range(K):
            gk_T = X[k] # 1 x N

            # 通过索引将gk_T,Ek收缩,不包含0
            I = [i for i,res in enumerate(gk_T) if res]

            # 全为0直接跳过不进行处理
            if len(I) == 0:
                continue

            #-- 求解收缩后的Ek
            D[:,k] = 0
            Ek_shrink = (Y - np.dot(D,X))[:,I]

            #-- 通过奇异值分解收缩后的Ek,得到新的dk和收缩后的gk_T
            #-- 奇异值分解,Sigma为奇异值向量,需要使用np.diag使其成为矩阵
            #-- 分解后的V需要进行转置再取第一列更新gk_T,Ek_shrink = U*Sigma*V_T
            U, Sigma, V = np.linalg.svd(Ek_shrink, full_matrices=False)

            #-- 用新的dk,gk_T 更新字典
            D[:,k] = U[:,0]
            X[k,I] = Sigma[0]*(V[0,:])
    return D, X

四、实验数据及数据预处理

1、试验数据预处理后保存为.npy格式

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
35
36
37
#-- 数据集目录
root_dir = r'DataSet'

#-- 用38张图片中的22张构建训练数据(8*8)
def save_train():
    files = os.listdir(root_dir)
    train_data = []
    width = 168
    height = 192
    for k in range(22):
        file = os.path.join(root_dir, files[k])
        img = plt.imread(file)
        for i in range(height//8):
            for j in range(width//8):
                train_data.append(img[i*8:i*8+8, j*8:j*8+8])

    #-- 将生成的训练数据用numpy保存
    train_data = np.array(train_data)
    np.save('train_data',train_data)

#-- 用38张图片中的后16张构建测试数据(8*8)
def save_test():
    files = os.listdir(root_dir)
    test_data = []
    width = 168
    height = 192
    for k in range(16):
        file = os.path.join(root_dir, files[k+22])
        img = plt.imread(file)
        for i in range(height//8):
            for j in range(width//8):
                test_data.append(img[i*8:i*8+8, j*8:j*8+8])

    #-- 将生成的训练数据用numpy保存
    test_data = np.array(test_data)
    np.save('test_data',test_data)

2、读取训练集和测试集并对测试集进行不同程度模糊(0,0.5,0.7)处理

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
#-- 读取22张训练集,number为读取的数量
def read_train(number:int=22):
    data = np.load('train_data.npy')
    return data[:number*504]

#-- 读取测试数据并进行不同程度模糊处理,number为读取的图片数量
#-- ratio为模糊个数
def read_test(ratio:int,number:int=16):
    if number > 22:
        print("测试数据不足!")
        return None

    data =  np.load('test_data.npy')
    data = data[:504*number]
    data = data.reshape(-1,64)

    # 每一行中随机选取 ratio个像素置0
    if not ratio:
        return data.reshape(-1,8,8)
    indexs = np.arange(64)
    for row in range(data.shape[0]):
        random_index = np.random.permutation(indexs)[:ratio]
        data[row][random_index] = 0
    return data.reshape(-1,8,8)

3、对预处理后的数据集进行逆向处理为原图片格式并显示

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
#-- 将处理后的数据转换为图像数据并显示
#-- test_data: -1 x 64
def reverse(test_data):
    test_data = test_data.reshape(-1,8,8)
    width = 168
    height = 192
    imgs = []
    for k in range(test_data.shape[0] // 504):
        img = np.zeros((height,width))
        for i in range(height//8):
            for j in range(width//8):
                img[i*8:i*8+8, j*8:j*8+8] = test_data[k*504 + i*21 + j]
        imgs.append(img)
    
    for i in range(len(imgs)):
        plt.subplot(5,5,i+1)
        plt.imshow(imgs[i], cmap='gray')
 
    plt.show()

# 显示原始训练集和测试集
def show():
    train = read_train().reshape(-1,64)
    reverse(train)
    test = read_test(0).reshape(-1,64)
    reverse(test)
    test = read_test(int(0.5*64)).reshape(-1,64)
    reverse(test)
    test = read_test(int(0.7*64)).reshape(-1,64)
    reverse(test)


效果
1、原始数据集

image-20191220105025486

2、原始测试集

image-20191220105041295

3、50%残差处理后的测试集

image-20191220105106148

4、70%残差处理后的测试集

image-20191220105125314

五、实验结果及实验分析

1、训练字典(迭代5次,稀疏度10)

#-- 对矩阵进行列归一化操作
#-- 使得列向量为单位向量
def normalization(A):
    for j in range(A.shape[1]):
        res_sum = np.sqrt(np.sum(np.square(A[:,j])))
        A[:,j] = np.divide(A[:,j],res_sum)
    return A

#-- 训练字典
#-- Y = 64 x 11088
#-- D = 64 x 441
#-- 稀疏编码 OMP
def train():
    S = 10  #  稀疏度
    L = 5   #  迭代次数
    e = 0.00001
    K = 441
    Y = read_train(10).reshape(-1,64).T

    #  对训练数据进行归一化处理
    Y_norm = np.zeros(Y.shape)

    #  每一列归一化
    for i in range(Y.shape[1]):
        normlize = np.linalg.norm(Y[:,i])
        mean = np.sum(Y[:,i])/Y.shape[0]
        Y_norm[:,i] = (Y[:,i]-mean)/normlize

    # -- 随机生成初始字典D0
    indexs = np.arange(Y.shape[1])[:K]
    D0 = Y_norm[:,indexs]
    
    # -- 归一化为单位向量
    D0 = normalization(D0)

    # 训练字典
    D, X = K_SVD(D0.copy(),S,Y_norm,e,L)

    #  保存结果
    np.save('D_train',D)
    np.save('X_train',X)

2、测试

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
35
36
37
38
39
40
41
42
43
44
#-- 利用学习的字典和稀疏矩阵对测试集进行重构
def test():
    D_file = 'D_train.npy'
    X_file = 'X_train.npy'

    D = np.load(D_file)
    X = np.load(X_file)

    # 获取模糊后的测试数据
    ratio = int(0.5 * 64)
    test_data = read_test(ratio) # -1 x 8 x 8
    test_data = test_data.reshape(-1,64).T # 64 x -1
    
    #  重构后的数据
    new_data = np.zeros(test_data.shape)

    # 归一化处理,只针对非零数据,然后使用OMP算法获取稀疏表达,计算后加入重构数据
    for i in range(test_data.shape[1]):
        index = np.nonzero(test_data[:,i])[0]
        if not index.shape[0]:
            continue
        normlize = np.linalg.norm(test_data[:,i][index])
        mean = np.sum(test_data[:,i][index])/index.shape[0]
        data_norm = (test_data[:,i][index]-mean)/normlize
        x = OMP(D[index,:],data_norm,0.0001,10)

        #  重构后需要反归一化处理
        new_data[:,i] = D.dot(x)*normlize + mean

    #  显示重构图像
    reverse(new_data.T)

    #  打印平均误差
    print(RMSE(read_test(0).reshape(-1,64).T,new_data))

#  计算原始数据和重构数据的误差
#  数据大小为 64 x -1
def RMSE(origin_data, constructed_data):
    error = 0
    for i in range(origin_data.shape[1]):
        error += np.sqrt(np.sum(np.square(origin_data[:,i]-constructed_data[:,i]))/64)
    
    return error / origin_data.shape[1] / 16

恢复效果
1、50%残差恢复效果

image-20191220105331171

2、70%残差恢复效果

image-20191220105355360

平均误差(RMSE)(5次迭代,稀疏度10)
1、50%残差RMSE: 0.2456584943753804
2、70%残差RMSE: 0.37693618384019295

六、小结

本次实验使用python复现了OMP算法和K-SVD算法,并完成了给定数据集的划分和预处理,通过训练字典,在测试集上重构残差图像,实现了人脸残差图像的恢复,并通过进一步的对比试验分析,了解了这些算法的一些性质(可能不完善不准确),基本达到了实验目的。

通过本课程的学习,了解了稀疏表达和低秩表达的相关概念以及非线性编码的一些概念,同时也了解到关于这些问题的一些解决方法。对于老师在课堂上的一些定理、结论的详细证明,也让人受益匪浅,感慨良多。而这些东西在人脸数据集上的实际应用也让人印象深刻,对表达学习这个领域也有了更多的好奇与兴趣。

参考文献

[1] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry and Y. Ma, “Robust Face Recognition via Sparse Representation,” in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 31, no. 2, pp. 210-227, Feb. 2009.

[2] 课程讲义

Read more »
叶藏也不行

叶藏也不行

叶藏也不行

4 posts
3 tags
RSS
github
Friends
  • jekyll
  • Milky Way
  • Lester
  • Herixth
  • Joke-Lin
© 2019 - 2024 叶藏也不行
Powered by Jekyll
Theme - NexT.Muse