网易笔试题：混合颜料


下面 

int getNumOfLeastColors(set<int>& colorSet) {// 求二进制矩阵的秩，即消元，最后看斜对角线上有几个 1  





的方法，就是求矩阵的秩


你就是一个画家！你现在想绘制一幅画，但是你现在没有足够颜色的颜料。为了让问题简单，我们用正整数表示不同颜色的颜料。你知道这幅画需要的n种颜色的颜料，你现在可以去商店购买一些颜料，但是商店不能保证能供应所有颜色的颜料，所以你需要自己混合一些颜料。混合两种不一样的颜色A和颜色B颜料可以产生(A XOR B)这种颜色的颜料(新产生的颜料也可以用作继续混合产生新的颜色,XOR表示异或操作)。本着勤俭节约的精神，你想购买更少的颜料就满足要求，所以兼职程序员的你需要编程来计算出最少需要购买几种颜色的颜料？ 
输入描述:

第一行为绘制这幅画需要的颜色种数n (1 ≤ n ≤ 50)

第二行为n个数xi(1 ≤ xi ≤ 1,000,000,000)，表示需要的各种颜料.

输出描述:

输出最少需要在商店购买的颜料颜色种数，注意可能购买的颜色不一定会使用在画中，只是为了产生新的颜色。

输入例子:

29

4096 8192 16384 32768 65536 131072 262144 524288 1048576 16 32 64 128 256 512 1024 2048 2097152 4194304 8388608 16777216 33554432 67108864 134217728 268435456 536870912 999999999 1000000000 15
输出例子:

27




思路：将所有的数二进制展开，构成一个矩阵，通过消元只保留对角线上 1，即求矩阵的秩，最后得到 1 的个数即秩的大小，也就是最少需要的颜色数。




C++实现：





[cpp] view
 plain copy





#include <iostream>  

#include <set>  

#include <vector>  

  

const int N = 8 * sizeof(int);// 整型 int 的二进制位数  

  

using namespace std;  

  

void getBinaryMat(const set<int>& s, vector<vector<int> >& bMat) {// 求每个数的二进制，存储在vector里面，最后得到二进制矩阵  

    if(s.empty()) return;  

    for(set<int>::iterator it=s.begin(); it != s.end(); it++) {  

        vector<int> v(N);  

        int flag = 1;  

        for(int i=0; i < N; i++) {  

            v[i] = (*it) & flag;  

            flag <<= 1;  

        }  

        bMat.push_back(v);  

    }  

}  

  

int getNumOfLeastColors(set<int>& colorSet) {// 求二进制矩阵的秩，即消元，最后看斜对角线上有几个 1  

    if(colorSet.empty()) return -1;  

    vector<vector<int> > bMat;  

    getBinaryMat(colorSet, bMat);  

    const int Row = bMat.size();  

    int r = 0, c = 0;  

    for(c=0; c < N && r < Row; c++, r++) {// 从每一列开始，将每一列消到只有 1 个 1 或者全 0  

        int i = 0;  

        for(i = r; i < Row; i++) {// 寻找这一列第一个非 0 的行  

            if(bMat[i][c] != 0)  

                break;  

        }  

        if(Row == i)  

            --r;  

        else {  

            swap(bMat[r], bMat[i]);  

            for(int k=i+1; k < Row; k++) {  

                if(0 != bMat[k][c]) {  

                    for(int j=c; j < N; j++) {  

                        bMat[k][j] ^= bMat[r][j];// 用第 r 行的 1 消除这一列上其他所有行的 1  

                    }  

                }  

            }  

        }  

    }  

    return r;  

}  

  

int main() {  

    int n;  

    while(cin >> n) {  

        set<int> colorSet;  

        int color;  

        for(int i=0; i < n; i++) {  

            cin >> color;  

            colorSet.insert(color);  

        }  

        cout << getNumOfLeastColors(colorSet) << endl;  

    }  

    return 0;  

}  

机器学习中会用到大量的数学操作，而 Numpy 计算库使这些操作变得简单，这其中就涉及到了 Numpy 的矩阵操作，下面我们就来一起学习如何在 Numpy 科学计算库中进行矩阵的一些基本运算。
1 矩阵的定义
定义矩阵使用 Numpy 科学计算库中的 mat 函数，如下所示：
numpy.mat(data, dtype=None)
data，表示矩阵的数据。
dtype，表示矩阵中的数据类型，默认是浮点数。
应用示例：
# (1) 定义一个3 X 3的矩阵，数据类型为 int
import numpy as np
data = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
A = np.mat(data, int)
print(A, type(A)) 

# 输出结果：
#[[1 2 3]
# [4 5 6]
# [7 8 9]] <class 'numpy.matrix'>

# (2) 定义一个3 X 3的矩阵，矩阵元素全为0，数据类型为 int
import numpy as np
A = np.mat(np.zeros((3, 3)), int)
print(A)

# 输出结果：
#[[0 0 0]
# [0 0 0]
# [0 0 0]]

# (3) 定义一个3 X 3的矩阵，矩阵元素全为1
import numpy as np
A = np.mat(np.ones((3, 3)))
print(A)

# 输出结果：
#[[1. 1. 1.]
# [1. 1. 1.]
# [1. 1. 1.]]

# (4) 定义一个3 X 3的单位矩阵
import numpy as np
A = np.mat(np.eye(3, 3), int)
print(A)

# 输出结果：
#[[1 0 0]
# [0 1 0]
# [0 0 1]]

# (5) 定义一个3 X 3的对角矩阵，主对角线之外的元素皆为0
import numpy as np
A = np.mat(np.diag([1, 2, 3]), int)
print(A)

# 输出结果：
#[[1 0 0]
# [0 2 0]
# [0 0 3]]

# (6) 定义一个3 X 3的矩阵，把100作为所有元素初始值
import numpy as np
A = np.mat(np.full((3, 3), 100), int)
print(A)

# 输出结果：
#[[100 100 100]
# [100 100 100]
# [100 100 100]]
2 矩阵的线性代数运算
2.1 矩阵的加法与减法
只有两个矩阵的行数和列数相等时，才可以进行矩阵的加法和减法运算，否则程序会抛出 ValueError 异常。
应用示例：
# 定义两个矩阵 A 和 B，分别进行矩阵的加法和减法运算
import numpy as np
A = np.mat(np.full((3, 3), 100), int)
B = np.mat(np.full((3, 3), 200), int)
print(A+B)
print(A-B)

# 输出结果：
#[[300 300 300]
# [300 300 300]
# [300 300 300]]
#[[-100 -100 -100]
# [-100 -100 -100]
# [-100 -100 -100]]
2.2 矩阵的数乘
某个实数乘以矩阵称作矩阵的数乘。
应用示例：
# 定义矩阵 A 和浮点数 a，进行矩阵的数乘运算
import numpy as np
a = 0.1
A = np.mat(np.full((3, 3), 100), int)
print(a*A)

# 输出结果：
#[[10. 10. 10.]
# [10. 10. 10.]
# [10. 10. 10.]]
2.3 矩阵的点乘
只有在第一个矩阵的列数与第二个矩阵的行数相等时，两个矩阵才能相乘，否则程序会抛出 ValueError 异常。
应用示例：
# 定义矩阵 A 和矩阵 B，进行矩阵的乘法运算
import numpy as np
A = np.mat(np.full((2, 3), 10), int)
B = np.mat(np.full((3, 3), 10), int)
print(A*B)  # 求矩阵相乘形式一
print(A.dot(B))  # 求矩阵相乘形式二
print(np.dot(A, B))  # 求矩阵相乘形式三

# 输出结果：
#[[300 300 300]
# [300 300 300]]
#[[300 300 300]
# [300 300 300]]
2.4 矩阵的转置
把矩阵的每一行转换为列，称为矩阵的转置。
应用示例：
# 定义矩阵 A，进行矩阵转置运算
import numpy as np
data = [[1,2,3],[4,5,6],[7,8,9]]
A = np.mat(data, int)
print(A.T)

# 输出结果：
#[[1 4 7]
# [2 5 8]
# [3 6 9]]
2.5 矩阵的求逆
非奇异矩阵下，可以对矩阵进行求逆运算。(非奇异矩阵就是行列式不为 0 的矩阵)
应用示例：
# 定义矩阵 A，进行矩阵求逆运算
import numpy as np
data = [[1, 2], [3, 4]]
A = np.mat(data, int)
print(A.I)

# 输出结果：
#[[-2.   1. ]
# [ 1.5 -0.5]]
2.6 矩阵的行列式
对于矩阵 A，均可对应一个标量 det(A)，它的值将告诉我们矩阵是否为非奇异的。
应用示例：
# 求矩阵 A 的行列式 det(A)
import numpy as np
A = np.mat([[1, 2],[3, 4]], int)
det = np.linalg.det(A)
print(det)

# 输出结果：
# -2.0000000000000004
计算结果不等于-2，是因为浮点数运算存在精度损失。
2.7 矩阵的秩
如果把矩阵看成一个向量组，那么秩就是线性无关向量的个数，也就是向量组的维度，概念比较复杂，有兴趣的读者可以继续探索。矩阵的秩应该是小于等于行数与列数的最小值。
# 求矩阵 A 的秩 rank(A)
import numpy as np
A = np.mat([[1, 2],[3, 4]], int)
rank = np.linalg.matrix_rank(A)
print(rank)

# 输出结果：
# 2
2.8 矩阵特征值和特征向量
A 为 n 阶矩阵，若数 λ 和 n 维非0列向量 x 满足 Ax=λx，那么数 λ 称为 A 的特征值，x 称为 A 的对应于特征值 λ的特征向量。
应用示例：
# 求矩阵 A 的特征值和其对应的特征向量
import numpy as np
A = np.mat([[1, 2],[3, 4]], int)
value, vector = np.linalg.eig(A)
print(value)
print(vector)

# 输出结果：
#[-0.37228132  5.37228132]
#[[-0.82456484 -0.41597356]
# [ 0.56576746 -0.90937671]]
2.9 矩阵的线性方程解
求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量。
应用举例：
# 求解如下线性方程组的解：
# x + y + z = 3
# 3x + y + 4z = 8
# 8x + 9y + 5z = 22
import numpy as np
A = np.mat([[1, 1, 1], [3, 1, 4], [8, 9, 5]], int)
b = np.mat([[3], [8], [22]], int)
x = np.linalg.solve(A,b)
print(x)

# 输出结果：
#[[1.]
# [1.]
# [1.]]

秩才是矩阵真正的“大小”。看看这个矩阵：
， 它的第二行和第三行都是第一行的倍数，熟知矩阵乘法行的观点我们就可以很容易把 
 表示为：
 。 
 的简化行阶梯矩阵 
 ，通过 
 我们很容易发现 
 的秩等于 
 。
秩 
 矩阵是很值得研究的，未来的学习中我们会发现它的行列式和特征值都很有意思。现在我们先来关注以下三点：
秩 
 矩阵只有一个主元。秩 
 矩阵只有一个主元列，也只有一个主元行，矩阵每列/行都是第一列/行的倍数。对于矩阵，列的性质总能决定一些行的性质，这是了不起的事实。
秩 
 矩阵可以表示为 
 ，反之能够表示为 
 的矩阵也是秩 
 矩阵。
注意：列矢量乘行矢量的结果是矩阵；点乘是行矢量乘列矢量，点乘的结果是一个数。另：我们讨论的矢量都是列矢量，如本例中的 
 ；当我们需要列矢量乘行矢量时，我们要把列矢量转置： 
 。
所有矩阵都可以表示成秩 
 矩阵的和。秩为 
 的矩阵可以表示为 
 个秩 
 矩阵的和。
接下来我们讲讲上面提到的第三点：秩为 
 的矩阵可以表示为 
 个秩 
 矩阵的和。我们借助消元过程实现这一目的。
将 
 简化成 
 的过程可以表示为： 
 ，所有消元的步骤相乘的结果是一个可逆矩阵 
，我们在第二章中仔细讲过这个结论。这个 
 我们可以通过 
 得到，将 
 变为 
 的行变换同时将 
 变成了 
 。这是我们在第二章中学过的Gauss-Jordan求矩阵的逆，不同点在于第二章中我们研究的是可逆方阵 
 ，它的简化行阶梯矩阵 
 为 
 ，因此消元矩阵 
 就等于 
 ，即 
 ，也就是 
 ，而这里的 
 可以是任意矩阵。
随后，我们自然而然的得到 
 ，下面我们看看这个矩阵分解是如何提供了一种用 
 个秩 
 矩阵和表示矩阵的方法的。我们用具体的例子来说明这个方法： 
 ，它的简化行阶梯矩阵 
 。 我们把 
 分解为
 。为什么最后这个等号能成立呢？因为根据矩阵相乘的第四种观点，也就是外积的观点， 
 ， 
 中的全零行自然是无足轻重的。实际上这里已经可以看到 
 被表示成了两个秩 
 矩阵的和了： 
 ，每个外积 
 都是一个秩 
 矩阵，而外积的总数由 
 非零行的数目决定，也就是矩阵的秩，这个例子中 
 的秩为
 。我们也可以写的再清晰一点 ：
 。
最后，值得一提的是：
 中 
 矩阵保留了 
 的主元列
 。这是由 
 的形式决定的， 
 中隐藏了 
 ，它选择了 
 的主元列。

一曲新词酒一杯，去年天气旧亭台。夕阳西下几时回？ 无可奈何花落去，似曾相识燕归来。小园香径独徘徊。 ———《浣溪沙·一曲新词酒一杯》——晏殊
更多精彩内容请关注微信公众号 “优化与算法”
 上一期介绍了低秩矩阵填充问题，这一期介绍一下低秩稀疏矩阵恢复问题。
1. 低秩矩阵恢复
将一个矩阵 
 分解为一个低秩矩阵部分 
 和一个独立同分布的高斯矩阵 
 的问题是经典的主成分分析（PCA）问题，可以通过对矩阵 
 进行奇异值分解得到最优解。
然而，当矩阵 
 为稀疏的噪声矩阵时，PCA不再适用于解决这个问题。此时 ，将一个矩阵 
 分解为一个低秩矩阵部分 
 和一个稀疏矩阵部分 
 的问题可以建模为下述优化问题：
其中 
，
 是观测矩阵。(1)式中 
 和 
 都是非线性且非凸的，优化起来非常困难，这个问题也被称为主成分追踪（Principal component pursuit, PCP）。幸运的是我们提前知道一些先验信息，即矩阵 
 是低秩的且矩阵 
 是稀疏的，从上一期介绍的关于矩阵填充理论中可知，矩阵的秩和 
 范数问题都可以进行凸松弛，从而为求解上述问题提供了途径。由于矩阵的核范数是矩阵秩的凸包络，矩阵的（1,1）范数是矩阵0范数的凸包络，因此可以将问题(1)松弛为如下凸优化问题：
求解式(2)也称为鲁棒主成分分析(RPCA)。
文献[1]中指出，只要低秩矩阵 
 的奇异值分布合理且稀疏矩阵的非零元素均匀分布，那么凸优化问题PCP就能够以接近1的概率从未知的任意误差中恢复出原始低秩矩阵 
 来。
求解(2)式的算法可以分为如下几大类：
迭代阈值算法 对于PCP问题时，迭代阈值算法(Iterative Thresholding, IT) 通过交替更新矩阵 
 和 
 来求解。IT算法的迭代形式简单且收敛，但它的收敛速度比较慢，且难以选取合适的步长。因此，IT算法具有有限的应用范围。
加速近端梯度算法 加速近端梯度算法(Accelerated Proximal Gradient, APG)的主要思想是利用了Nesterov加速，使算法能够快速收敛。
对偶方法 将原问题转化为其对偶问题(非线性、非光滑)，并使用最速上升法等可以求解。对偶方法比APG算法具有更好的可扩展性，这是因为在每次迭代过程中对偶方法不需要矩阵的完全奇异值分解。
增广拉格朗日乘子法
这些方法都非常经典，这里不再细述，总的来说，只要将问题转化为凸问题，就有一大堆方法可以用来求解。这里仅介绍一种增广拉格朗日乘子算法，即交替方向方法(Alternating direction methods, ADM)，也称为不精确拉格朗日乘子法(Inexact ALM, IALM)。 下面给出上述几种算法的比较(数据来源于网络) 
2. 交替方向算法(ADM)
对于优化问题(2)，首先构造增广拉格朗日函数：
当 
 时，使用交替方法求解块优化问题：
使用精确拉格朗日乘子法(Exact ALM, EALM)交替迭代矩阵 
 和 
，直到满足终止条件为止。若 
，则
再根据 
 更新矩阵 
：
记 
 和 
 分别为 
 和 
 的精确值，则矩阵 
 的更新公式为：
参数 
 可以更新如下：
其中 
 为常数，
 为一个小的正数。
上述精确ALM方法在内循环中要多次更新，进行多次奇异值分解，为此文献[1]提出了非精确拉格朗日乘子法(Inecact ALM, IALM)，它不需要在每次外循环开始之前要求 
 的精确解，也就是去掉了ALM方法的内循环，其更新公式变成了如下形式：
上面式子中的奇异值阈值算子 
 和软阈值算子 
 的定义参见上一期<低秩矩阵填充|奇异值阈值算法>。
##4. 低秩矩阵恢复的应用 低秩矩阵恢复技术是一个非常有研究价值和实用价值的技术，它的应用也非常广泛，比如说：
视频背景建模。 
图像恢复(去光照、阴影等) 
图像类别标签净化
文本主题分析
音乐词曲分离
图像矫正与去噪 
图像对齐
##5. 仿真
ADM算法matlab代码如下：
function [L,S] = pcp_ad(M,u,lambda,itemax,tol)
% solve PCP problem by ADM algorithm
% v1.0 2020-1-1
% function：solve the following optimization problem
%                  min  ||X||*+lambda||E||_F
%                  s.t. M = A+E

% initialize S0 and Y0 and L0
[m,n] = size(M) ;
S = zeros(m,n) ;
Y = S ;
L = S ;

% the observed matrix can contain non number
unobserved = isnan(M);
M(unobserved) = 0;

if nargin < 2
    lambda = 1 / sqrt(max(m,n));
end
if nargin < 3
    u = 10*lambda;
end
if nargin < 4
    tol = 1e-6;
end
if nargin < 5
    itemax = 1000;
end

for ii = 1:itemax
    L = sig_thre(M-S+(1/u)*Y,(1/u)) ;
    S = soft_thre(M-L+(1/u)*Y, lambda/u) ;
    Z = M-L-S ;
    Y = Y+u*Z ;
    error = norm(M-L-S,'fro')/norm(M,'fro') ;
    if (ii == 1) || (mod(ii, 10) == 0) || (error < tol)
        fprintf(1, 'iter: %04dterr: %ftrank(L): %dtcard(S): %dn', ...
            ii, error, rank(L), nnz(S));
    end
    if error<tol
        break;
    end
end

数值测试代码：

% solve PCP problem by alternating direction method
clear
clc

m = 100 ;
n = 100 ;
r = 0.05*n ;
rate = 0.05 ;
% Generating a low rank matrix
LL = randn(m,r)/sqrt(m)*randn(r,n)/sqrt(n) ;
% Generating a large sparse noise matrix (Bernoulli matrix)
SS = randi([0,1],m,n) ;
SS(SS==0) = -1 ;

% sampling
ss = SS(:) ;
index = sort(randperm(m*n,ceil(rate*n*m))) ;
ss1 = zeros(m*n,1) ;
ss1(index) = ss(index) ;
SS = reshape(ss1,m,n) ;
M = LL+SS ;

lambda = 1/sqrt(max(m,n)) ;
u = 10*lambda ;

% [L,S] = pcp_ad(M,u,lambda,1000) ;
[L,S] = RobustPCA(M,lambda,u);
% [L,S] = pcp_ad(M,u,lambda,500,1e-8);
% [L,S] = adm_lrr(M);
MM = M-L-S ;

norm(M-MM,'fro')/norm(M,'fro')
norm(M-L-S,'fro')/norm(M,'fro')
norm(L-LL,'fro')/norm(LL,'fro')
norm(S-SS,'fro')/norm(SS,'fro')

运行上面程序，显示结果norm(M-L-S,'fro')/norm(M,'fro')约为9e-7，norm(L-LL,'fro')/norm(LL,'fro')约为1e-5。
低秩图像恢复仿真程序：
% low rank and sparse noise image recovery
clear
clc

A = imread('C:xxxxxxxxx.bmp') ;

WW = double(A) ;
a1 = double(A(:,:,1)) ;
a2 = double(A(:,:,2)) ;
a3 = double(A(:,:,3)) ;
[M,N] = size(a1);
X = zeros(M,N,3);
    
for jj = 1:3
    lambda = 1*1 / sqrt(max(M,N)); 
    u =  1*lambda;
    [ X(:,:,jj),S(:,:,jj)] = RobustPCA(WW(:,:,jj),lambda,u,1e-8,200) ;
end

figure(1)
subplot(3,1,1)
imshow(A)
title("原图",'fontsize',12)
subplot(3,1,2)
imshow(uint8(X))
title("低秩图",'fontsize',12)
d = S ;
d(d<20) = 255 ;
subplot(3,1,3)
imshow(uint8(d))
title("噪声图",'fontsize',12)

低秩图像恢复结果如下： 
从上面图像恢复结果来看，效果还不错。
参考文献
[1] Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11. 
[2] 史加荣, 郑秀云, 魏宗田, & 杨威. (2013). 低秩矩阵恢复算法综述. 计算机应用研究, 30(6), 1601-1605. 
[3] Cui, X., Huang, J., Zhang, S., & Metaxas, D. N. (2012, October). Background subtraction using low rank and group sparsity constraints. In European Conference on Computer Vision (pp. 612-625). Springer, Berlin, Heidelberg. 
[4] Wright, J., Ganesh, A., Rao, S., Peng, Y., & Ma, Y. (2009). Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization. In Advances in neural information processing systems (pp. 2080-2088). 
[5] Peng, Y., Ganesh, A., Wright, J., Xu, W., & Ma, Y. (2012). RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE transactions on pattern analysis and machine intelligence, 34(11), 2233-2246.
更多精彩内容请关注微信公众 “优化与算法”